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In  view  of  the  increasing  interest  in  using  composite 
materials  for  aerospace  structures,  the  analysis  of 
laminated  composite  plates  becomes  essential.  A systematic 
study  of  the  response  of  laminated  plates  subject  to  static, 
dynamic , and  impact  loading  has  been  conducted  by  using  a 
hybrid  stress  finite  element  method. 

The  hybrid  stress  model  was  based  on  the  modified 
Reissner  principle.  The  transverse  shear  effects,  as  well 
as  the  coupling  effects  of  stretching  and  bending,  were 
included  in  the  model.  The  displacement  field  was 
interpolated  through  shape  functions  and  the  stress  field 
was  interpolated  through  assumed  stress  polynomials. 

The  validity  of  the  hybrid  stress  model  was  determined 
by  comparing  the  numerical  results  with  existing  exact 
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analytical  solutions.  A frontal  solution  method  was 
employed  to  solve  the  cylindrical  bending  problems.  A 
subspace  iteration  method  was  used  to  compute  the 
fundamental  frequency  of  a simply  supported  laminated  plate. 
The  Newmark  direct  integration  method  was  used  to  integrate 
the  finite  element  equations  of  motion  step  by  step. 
Excellent  accuracy  and  fast  convergence  were  observed  in  the 
numerical  results. 

The  dynamic  response  of  a laminated  plate  impacted  by  a 
projectile  has  been  analyzed.  The  force  distributions  were 
based  on  the  principle  of  linear  impulse  and  momentum  and 
the  Hertzian  impact  law.  The  transverse  deflection  at  the 
center  of  laminate  was  in  good  agreement  with  the  deflection 
determined  experimentally.  The  laminate  bending  stiffness, 
projectile  velocity,  and  laminate  thickness  had  effects  on 
the  deflection  and  stress  distributions.  The  deflection  and 
stress  distributions  showed  evidence  of  flexural  wave  motion 
induced  by  impact  on  the  laminated  plate.  The  flexural  wave 
motion  may  be  the  cause  of  experimentally  observed 
delamination  crack  propagation. 

The  interlaminar  transverse  shear  stresses  and  normal 
stress  were  incorporated  into  a quadratic  failure  criterion 
to  perform  the  failure  analysis.  The  results  showed  that 
the  delamination  crack  was  initiated  near  the  area  of  impact 
and  propagated  outwardly.  They  also  showed  that  the 
delamination  crack  for  both  interfaces  possibly  occurred 
simultaneously . 
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CHAPTER  1 


INTRODUCTION 

1.1  Motivation  for  the  Research 

During  the  past  two  decades,  the  applications  of 
composite  materials  have  grown  rapidly  and  are  now  very 
commonplace  around  the  world.  Industries  such  as  aircraft, 
automobiles,  sporting  goods,  electronics,  and  appliances  are 
quite  dependent  on  fiber-reinforced  composites.  In  the 
contemporary  technology  of  structural  composite  materials, 
major  deficiencies  still  exist  with  respect  to  the  ability 
to  determine  the  stress  field  within  a multilayered 
composite  laminate  [1],  especially  for  the  dynamic  response. 

All  known  approximate  laminate  theories  are  based  upon 
assumed  displacement  fields.  These  theories  lead  to  results 
which  cannot  describe  the  actual  response  of  laminates  [2, 
3,  4,5,  6,  7].  To  get  an  acceptable  laminate  field  theory 

some  basic  requirements  need  to  be  satisfied.  These 
requirements  are  stated  below. 

a)  All  six  stress  components  are  non-zero  in  general; 

b)  The  stresses  within  each  layer  satisfy  the 
equilibrium  equations; 

c)  Traction  and  displacement  continuity  conditions  at 
the  interface  between  adjacent  layers  are  satisfied. 
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The  hybrid  stress  model  proposed  here  has  fulfilled  the 
above  requirements.  For  this  reason,  the  author  believes 
that  the  hybrid  stress  finite  element  method  will  be  a more 
proper  approach  than  the  conventional  displacement  finite 
element  method  for  the  stress  analysis  of  laminated 
composite  plates. 

Impact  damage  in  laminated  composite  plates  has  been  a 
major  concern  in  engineering  applications.  Delamination 
failures  are  often  exhibited  at  the  interfaces  of  laminated 
plates  which  are  subject  to  local  impact.  The  understanding 
of  when,  where,  and  how  the  delamination  occurred  will  be 
very  helpful  to  the  engineering  design  of  composite 
structure  components. 

The  objective  of  this  research  was  to  use  the  hybrid 
stress  finite  element  method  to  investigate  the  stress  and 
displacement  fields  of  laminated  composite  plates  subject  to 
static,  dynamic,  and  impact  loading  and  the  natural 
vibration  frequencies  of  laminated  composite  plates.  A 
suitable  failure  criterion  is  proposed  and  incorporated  to 
examine  the  delamination  mechanisms  of  laminated  plates 
under  impact  loading.  The  results  of  this  study  will  lead  to 
more  concrete  concepts  about  delamination  mechanisms,  and 
this  study  is  expected  to  be  a starting  point  for  further 
understanding  of  delamination  mechanisms. 
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1.2  Laminated  Composite  Structure 
The  word  "composite"  means  "consisting  two  or  more 
distinct  parts."  A material  having  two  or  more  distinct 
constituent  materials  or  phases  is  considered  as  composite 
material.  Composite  materials  have  a long  history  of  usage. 
Their  beginnings  are  unknown,  but  all  recorded  human  history 
contains  references  to  some  form  of  composite  materials.  For 
example,  straw  was  used  by  the  Israelites  to  strengthen  mud 
bricks.  Plywood  was  used  by  the  ancient  Egyptians  when  they 
realized  that  the  wood  could  be  arranged  to  achieve 
superior  strength  and  resistance  to  thermal  expansion  [8]. 
More  recently,  fiber-reinforced  resin  composites,  which  have 
high  s trength-to-weight  and  stiff ness-to-weight  ratios,  have 
become  important  in  weigh t- sens i ti v e applications  like 
aircraft  and  space  vehicles. 

Material  properties  that  can  be  improved  by  forming  a 
composite  material  include  [8] 

. Strength  . Temperature-dependent  behavior 

. Stiffness  . Thermal  insulation 

. Corrosion  resistance  . Thermal  conductivity 

. Wear  resistance  . Acoustical  insulation 

. Weight  . Attractiveness 

. Fatigue  life  . Damping 

A lamina  usually  consists  of  a single  layer  of 
unidirectional  reinforcement  in  a matrix.  Several  laminae 
are  stacked  together  to  form  a structure  termed  a laminate 
as  shown  in  Figure  1.1.  The  principal  material  directions  of 
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(a)  Lamina 


Figure  1.1.  Lamina  and  laminate  construction 
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each  lamina  in  a laminate  have  various  orientations  and  the 
stacking  sequence  is  not  necessarily  symmetric  about  the 
midsurface  of  the  laminate.  The  layers  of  a laminate  are 
usually  bonded  together  by  the  same  matrix  material  used  in 
the  lamina.  The  bond  between  two  laminae  in  a laminate  is 
assumed  to  be  perfect;  thus,  the  laminae  cannot  slip  over 
each  other.  The  displacements  of  laminae  remain  continuous 
across  the  bond.  One  of  the  most  important  advantages  of 
fibrous  composite  laminates  is  that  their  anisotropic 
properties  can  be  controlled  very  effectively;  that  is, 
altering  the  fiber  directions  of  each  lamina  can  give  the 
laminate  different  strengths  and  stiffnesses  in  various 
directions.  Thus,  the  strengths  and  stiffnesses  of  a 
laminated  fiber-reinforced  composite  can  be  tailored  to 
match  the  specific  design  requirements  of  a structural 
element  being  built. 

A potential  problem  in  the  laminate  construction  is  the 
introduction  of  shearing  stress  at  the  interface  between  two 
laminae.  The  shearing  stresses  arise  due  to  the  tendency  of 
each  lamina  to  deform  independently  of  its  neighbors  because 
each  lamina  may  have  different  properties.  Such  interlaminar 
shear  stresses  as  well  as  normal  stress  are  the  major  causes 
for  the  delamination  mechanisms,  and  this  subject  has  been 
investigated  by  many  authors  in  recent  years. 
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1.3  Stress  Analysis  by  Using  the  Finite  Element  Method 

In  the  earliest  developments  of  finite  element  methods 
almost  all  emphasis  was  directed  toward  the  development  of 
effective  finite  elements  for  the  solution  of  specific 
structural  analysis  problems.  During  the  past  two  decades, 
the  broad  potential  of  the  method  was  rapidly  recognized  and 
more  general  techniques  for  structural  analysis  were 
developed  while  at  the  same  time  the  techniques  were  well 
established  within  many  scientific  and  engineering 
disciplines.  Although  there  may  be  considerable  diversity  in 
the  formulation,  a finite  element  method  can  be 

distinguished  by  the  following  features  [9]: 

(1)  The  physical  region  of  the  problem  is  subdivided 
into  subregions  or  finite  elements. 

(2)  One  or  more  of  the  dependent  variables  is 

approximated  in  functional  form  over  each  element 
and  hence  over  the  whole  domain.  The  parameters  of 
this  approximation  subsequently  become  the  unknowns 
of  the  problem. 

(3)  Substitution  of  the  approximations  into  the 
governing  equations  (or  their  equivalent)  yields  a 
set  of  algebraic  equations  in  the  unknown 
parameters.  The  solution  of  these  equations  yields 
the  parameters  and  hence  the  approximate  solution 
to  the  problem. 

Most  commonly,  the  governing  equations  are  not  used 
directly,  but  an  alternative  formulation  such  as  a 
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variational  principle  is  used  instead.  Since  the  number  of 
unknowns  in  the  final  set  of  equations  is  often  very  large, 
it  is  a common  practice  to  use  matrix  notation,  for  both 
conciseness  and  facilitation  of  computer  programming. 

In  1964,  the  assumed  hybrid  stress  model  was  suggested 
by  Pian  [10]  to  derive  the  element  stiffness  matrices  and 
was  later  extensively  applied  to  the  analysis  of  plate  and 
shells  as  well  as  crack  and  other  problems  by  many  authors. 
The  derivation  is  based  on  the  principle  of  minimum  comple- 
mentary energy.  In  this  formulation  the  displacements  of  the 
element  boundary  are  assumed  so  as  to  provide  displacement 
compatibility  between  elements,  and  internal  stress  varia- 
tions are  assumed  to  satisfy  the  differential  equations  of 
equilibrium.  Greater  flexibility  is  achieved  in  the  stress 
distribution  by  the  introduction  of  elemental  assumed  stress 
polynomials  [11]. 

1.4  Impact  on  a Laminated  Plate 

The  transient  response  of  laminated  plates  subject  to 
local  impact  loadings  has  been  of  significant  concern  in 
many  advanced  engineering  structures  and  components.  For 
example,  the  leading  edge  of  an  aircraft  wing  or  blades  in  a 
jet  engine  may  be  hit  by  foreign  objects  such  as  stones, 
nuts,  bolts,  hailstones  or  birds.  In  these  cases  the  impact 
resistance  of  laminated  plates  must  be  known  in  order  to 
design  the  laminated  plates  which  meet  aircraft  safety 
requirements . 
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Impact  problems  with  energies  far  below  the  penetration 
levels  can  cause  significant  damage  to  components  made  of 
composite  materials,  since  delamination  often  occurs  at 
energy  levels  which  are  an  order  of  magnitude  below  those 
causing  penetration  [12]. 

Plate  impact  experiments  have  been  conducted  at  the 
University  of  Florida.  In  the  experiments,  a small  steel 
cylinder  is  fired  from  a gas  gun  as  shown  in  Figure  1.2  at 
subperforation  velocities  of  10-80  m/second.  The  target  is  a 
fiber-reinforced  composite  plate  clamped  at  the  four  edges 
as  shown  in  Figure  1.3.  At  sufficiently  high  subperforation 
speeds  the  delamination  failure  is  usually  observed  at  the 
interface  of  two  adjacent  layers  of  different  orientation. 

In  the  present  study  a dynamic  hybrid  stress  finite 
element  program  is  developed  to  calculate  the  displacement 
and  stress  distributions  in  a composite  laminated  plate 
subject  to  a local  impact  loading.  Transient  response  of 
displacements  and  stresses  is  presented.  Interlaminar  shear 
stresses  and  transverse  normal  stress,  which  are  responsible 
for  delamination,  are  incorporated  into  a failure  criterion 
to  examine  the  delamination  crack  initiation  and 
propagation.  Although  the  results  predicted  herein  are  not 
expected  to  match  those  of  the  experimental  results  exactly, 
the  basic  characteristics  of  the  whole  model  have  the 
potential  to  describe  why  the  delamination  crack  occurs  and 
are  consistent  with  experimental  observations. 
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Figure  1.2.  Schematic  dravzing  of  gas  gun  assembly  and  related  equipments 
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Figure  1.3.  Plate  impacted  by  a projectile 
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1 . 5 Review  of  the  Literature 

Classical  Lamination  Theory  (CLT)  [2,  3]  has  wide  usage 
as  the  basis  for  the  stress  analysis  and  design  of 
structural  laminates.  Basic  assumptions  of  CLT  are  a)  The 
Kirchhoff  assumption  of  negligible  transverse  shear  strain 
and  transverse  normal  strain  constitutes  a statement  of 
undeformable  normals  to  the  midsurface  of  the  plate;  b)  The 
assumption  of  linear  in-plane  displacements  through  the 
thickness  and  the  transverse  displacements  are  each  constant 
through  thickness;  c)  The  assumption  that  in-plane  strains 
are  small  compared  to  unity  (small  strain  theory);  d)  Rotary 
inertia  effects  are  neglected.  Because  of  these 
characteristics,  which  limit  its  generality  in  the 
description  of  laminate  response,  CLT  can  predict  accurate 
stresses  and  deflections  only  when  plate  thickness  is  very 
small  while  compared  to  its  length  and  width  (thin  plate) 
and  the  deflection  is  small  compared  to  the  thickness. 

A modified  theory  for  isotropic  homogeneous  plates  was 
introduced  by  Mindlin  [13].  The  theory  includes  transverse 
shear  and  rotary  inertia  effects.  Later,  Yang,  Norris,  and 
Stavsky  [14]  extended  Mindlin 's  theory  to  analyze  laminated 
plates.  Because  of  a lack  of  solutions  to  specific  boundary- 
value  problems,  its  consequences  had  not  been  explored  until 
Whitney  and  Pagano  [4]  presented  a detailed  investigation  of 
the  theory  developed  in  [14].  They  solved  a number  of 
specific  boundary-value  problems  of  simply  supported 
laminated  plates.  Those  problems  were  related  both  to  static 


12 


bending  under  sinusoidal  loading  and  to  natural  vibration 
frequencies.  For  the  cylindrical  bending  problem,  they  found 
the  effects  of  shear  deformation  on  the  maximum  deflection 
are  very  significant  in  composites  with  a span-to-depth 
ratio  as  high  as  20.  The  stresses  predicted  by  this  theory 
were  identical  to  those  of  Classical  Lamination  Theory. 

Pagano  [15,  16]  and  Pagano  and  Hatfield  [17] 
constructed  3-dimensional  elasticity  solutions  of  a number 
of  specific  boundary  value  problems  under  cylindrical 
bending.  The  results  were  compared  with  the  CLT  solutions 
for  the  stresses  and  displacements.  It  revealed  that  CLT 
leads  to  a very  poor  description  of  laminate  response  at  low 
span-to-depth  ratios.  Pipes  and  Pagano  [18]  applied  the 
elasticity  theory  and  the  finite  difference  technique  to 
investigate  the  composite  laminate  response  under  uniform 
axial  strain.  The  results  of  the  application  showed  that 
significant  interlaminar  shear  stresses  are  required  to 
allow  shear  transfer  between  the  laminae  of  the  laminate. 
Applications  of  these  analytical  techniques  are  restricted 
to  some  special  cases.  Some  of  the  examples  are  plies  with 
special  orientations,  simply  supported  boundary  conditions, 
and  sinusoidally  distributed  loading.  In  practice,  more 
complicated  structural  members  are  presented  with  holes, 
curved  shapes,  edge  loading,  and  discontinuities  in 
thickness.  These  complex  factors  have  limited  the  attempts 
to  use  analytical  techniques,  which  leads  to  the  dependence 
on  numerical  methods  and  digital  computers. 
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Finite  element  methods  have  been  widely  applied  in 
composite  structure  analyses.  Pryor  and  Barker  [5]  suggested 
a laminated  plate  element  including  transverse  shear  effects 
to  analyze  the  static  response  of  a laminated  plate  in 
cylindrical  bending.  The  theory  was  based  on  the  assumption 
that  the  deformed  normal  remains  straight  but  is  not  normal 
to  the  midsurface.  Poor  descriptions  for  the  response  of  a 
three-layer  cross-ply  thick  laminate  were  noticed.  The 
predicted  results  were  very  close  to  the  CLT  results.  Putcha 
and  Reddy  [19]  and  Reddy  and  Chao  [20]  developed  a shear 
flexible  finite  element  based  on  the  Hencky-Mindlin  shear 
deformation  theory.  The  behavior  of  a laminated  plate  in 
static  bending  was  investigated.  Some  disagreements  were 
found  between  the  mixed  shear  flexible  finite  element 
solution  and  3-dimensional  elasticity  solution.  The 
difference  was  more  severe  in  a thick  plate.  A higher  order 
shear  deformation  theory  of  a laminated  composite  plate, 
based  on  the  shear  deformation  theory  of  Whitney  and  Pagano 
[4],  but  accounting  for  a parabolic  distribution  of  the 
transverse  shear  strains  through  the  thickness  of  the  plate 
was  developed  by  Reddy  [21].  This  method  predicted  the 
deflections  and  stresses  more  accurately  than  [20],  but  some 
differences  still  existed  in  the  thick  plate  results  when 
compared  with  Pagano 's  elasticity  results  [16]. 

Since  Pian  [10]  first  established  the  assumed  stress 
hybrid  finite  element  model  and  derived  the  corresponding 
element  stiffness  matrix  in  1964,  he  and  his  colleagues  have 
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shown  that  such  a hybrid  model  not  only  provides  the 
flexibility  and  ease  in  fulfilling  the  compatibility 
condition  of  the  finite  element  method  for  plates  and  shells 
but  is  also  highly  accurate.  In  the  meantime,  many  other 
authors  have  also  confirmed  the  usefulness  of  such  a model 
and  have  adopted  it  for  many  applications  in  solid  mechanics 
as  well  as  fluid  mechanics.  The  following  is  a survey  of 
some  of  the  works  that  have  been  investigated  by  using  the 
hybrid  stress  finite  element  method.  Pian  and  Tong  [22]  used 
a variational  principle  to  derive  the  hybrid  model.  They 
solved  plate  bending  problems  including  both  thermal  effects 
and  transverse  shear  effects.  Laminated  thick  plate  static 
bending  problems  were  solved  by  Mau , Tong,  and  Pian  [23].  In 
the  comparison  of  results  with  the  elasticity  solution  [15, 
16],  they  observed  excellent  accuracy  in  predicting  both 
deflections  and  stresses.  In  their  assumption  for  the  stress 
field,  transverse  normal  stress  was  not  included.  Constant 
transverse  displacement  through  the  thickness  of  the 
laminate  was  also  assumed.  These  facts  did  not  agree  with 
the  actual  mechanism  of  deformation  of  laminated  plate 
bending.  Tong,  Pian,  and  Lasry  [24]  used  a hybrid  model  to 
solve  crack  problems  in  plane  elasticity. 

Spilker  [25]  and  Spilker  and  Munir  [26]  formulated  a 
hybrid  stress  element  model  for  thick  multilayer  laminate 
with  the  assumption  of  higher  order  through-the-thickness 
distributions  for  stress  and  displacement  within  each  ply. 
Very  good  correlations  between  computed  and  elasticity 
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solutions  were  shown,  but  it  was  restricted  to  cylindrical 
bending  of  cross-ply  laminate  and  two-dimensional  plane 
strain  element.  The  analysis  of  axisymmetric  structures 
subject  to  non-axisymmetric  loading  was  investigated  by 
Spilker  and  Daugirda  [27].  In  a series  of  papers,  Spilker 
and  his  associates  [28,  29,  30,  31]  derived  alternative 
three-dimensional  and  plane  isoparametric  hybrid  stress 
elements  for  thick  and  thin  plates;  each  had  different 
assumptions  for  stress  and  displacement  fields.  These  models 
were  not  very  general  since  there  were  usually  some 
limitations  imposed.  Later,  Spilker  [32,  33]  developed  an 
eight-node  isoparametric  multilayered  plate  element  for  the 
analysis  of  thin  to  thick  fiber-reinforced  composite  plates. 
The  new  model  was  more  general  than  previous  models  in 
describing  the  laminated  plate  response,  but  still  retained 
the  assumption  of  constant  transverse  displacement  through 
the  laminate  thickness.  Some  static  bending  problems  and 
computation  cost  for  this  model  were  discussed  in  these 
papers . 

Rhee  and  Atluri  [34]  presented  the  solution  of  bending 
stress  intensity  factors  of  a plate  with  a through-the- 
thickness  crack.  Analyses  of  large  deformation  of  inelastic 
solids  based  on  hybrid  stress  finite  elements  were 
undertaken  by  Reed  and  Atluri  [35].  Horrigmoe  [36]  described 
a hybrid  stress  finite  element  formulation  for  geometrically 
non-linear  analysis  of  thin  shell  structures.  Tong  and  Fung 
[37]  investigated  the  slow  particulate  flow  in  channels  and 
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tubes.  Rubinstein  et  al . [38]  have  discussed  the  hybrid 
stress  elements  in  a theoretical  way.  A method  of  selection 
of  a least-order,  stable,  invariant  stress  field  was 
proposed . 

Free  vibrations  and  dynamic  response  of  laminated 
structures  were  investigated  by  many  authors.  Noor  [39] 
presented  the  fundamental  frequencies  of  multilayered 
composite  plates  based  on  3-dimensional  elasticity  theory. 
Pillasch,  Majerus,  and  Zak  [40]  developed  a finite  element 
model  for  the  dynamic  analysis  of  laminated  thick  plates.  A 
three-layer  simply  supported  laminated  plate  subject  to  an 
impulsive  loading  was  examined.  Numerical  results  of  the 
laminate  transverse  deflections  were  presented.  By  using  the 
shear  flexible  finite  element  method,  Reddy  [41,  42,  43,  44] 
and  Putcha  and  Reddy  [45]  made  a very  detailed  investigation 
on  the  free  vibration  and  transient  response  of  isotropic, 
orthotropic  and  layered  anisotropic  composite  plates. 
Different  types  of  loadings , boundary  conditions  , and 
orientation  of  each  lamina  were  considered  in  his  sample 
problems . 

Mau,  Pian,  and  Tong  [46]  used  a hybrid  stress  element 
with  an  appropriately  computed  mass  matrix  to  investigate 
the  vibration  of  laminated  plates  and  shells.  In  the  same 
year,  they  [47]  suggested  a unifying  approach  in  deriving 
the  geometric  stiffness  and  mass  matrices  for  hybrid  finite 
element  models.  The  variational  formulation  was  based  on  a 
modified  Reissner  principle.  A simple  beam  example  was 
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presented  in  this  paper.  Spilker  et  al . [48]  adopted  Mau  and 
others'  [23]  quadrilateral  element  and  stress  field 
assumption  to  investigate  static  and  dynamic  bending 
problems.  Wang  [49]  employed  Tong  and  others'  [47]  approach 
for  the  derivation  of  mass  and  stiffness  matrices  to  develop 
a finite  element  method  together  with  hybrid  stress  crack- 
tip  elements  to  analyze  the  delamination  mechanism  of  a 
laminated  plate  under  dynamic  loading. 

Concerning  the  impact  problems,  Goldsmith  [50]  provided 
detailed  discussion  for  impact  on  isotropic  beams  and 
plates.  Chow  [51]  investigated  the  transient  response  of  a 
rectangular  laminated  plate  subject  to  a concentrated  impact 
loading  by  using  shear  deformation  plate  theory  [4].  Chen 

[52]  applied  a Hertzian  contact  law  to  analyze  the  response 
of  impact  on  anisotropic  materials.  Sun  and  Chattopadhyay 

[53]  developed  a nonlinear  integral  equation  to  solve  for 
the  contact  force  and  the  dynamic  response  of  a simply 
supported  laminate  subject  to  central  impact.  The  dynamic 
response  of  laminates  impacted  by  spherical  impactors  was 
discussed  by  Greszczuk  and  Chao  [54]  using  Hertzian  contact 
theory.  The  delamination  mechanisms  in  a laminated  composite 
plate  produced  by  the  impact  of  a blunt-ended  impactor  were 
investigated  by  Cristescu,  Malvern,  and  Sierakowski  [55].  A 
systematic  experimental  study  of  failure  mechanisms  on 
impact  of  laminated  plates  was  conducted  by  Takeda  [56].  The 
effect  of  impactor  mass  and  velocity,  and  ply  orientation  on 
laminated  plate  delamination  was  investigated  by  Takeda, 
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Sierakowski , and  Malvern  [57].  Sun,  Sankar , and  Tan  [58] 
established  a contact  law  for  a composite  plate  and  a steel 
ball.  The  contact  behavior  was  incorporated  into  a finite 
element  program  to  analyze  the  wave  propagation  phenomenon. 
Moyer  and  Ghashghai-Abdi  [12]  employed  the  finite  element 
method  and  finite  difference  method  to  investigate  the 
dynamic  response  of  beam  and  plate  components  subject  to 
impact  loading  with  geometrically  nonlinear  deformation. 
Keer  and  Woo  [59]  used  elasticity  theory  and  plate  theory  to 
determine  the  contact  force  history  and  pressure 
distribution  on  a circular  plate  with  low  velocity  impact. 


CHAPTER  2 


FORMULATION  OF  THE  HYBRID  STRESS 
FINITE  ELEMENT  METHOD 

2.1  Definition  of  Material  Elastic  Constants 
The  strain-stress  relations  for  anisotropic  materials 
can  be  expressed  as 

e i = sijaj  if  j = If  6 2.1 

where  are  the  strain  components,  Sj_j  are  the  elements  of 
the  compliance  matrix,  and  are  the  stress  components.  For 
a unidirectional ly  fiber-reinforced  lamina  with  axes  x3  , 
x2 , x3  aligned  with  the  principal  material  directions  as 
shown  in  Figure  2.1  (x-j_  corresponds  to  fiber  direction,  x2 
is  the  in-plane  coordinate  perpendicular  to  x^  , and  x3  is 
the  transverse  coordinate  which  is  perpendicular  to  the  x^- 
x 2 plane),  the  strain-stress  relations  for  an  orthotropic 

lamina  are 
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The  six  stresses  are  shown  in  Figure  2.2  etc.,  where  T23=a4 ' 
x13=a5'  t12=  a6’  The  comPliance  components  in  terms  ofthe 
engineering  constants  are  described  as 
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E^,  E2  , E3  = Young's  moduli  in  1 , 2,  and  3 directions, 
respectively . 

= Poisson's  ratio  for  transverse  strain  in 
the  Xj  direction  when  stressed  in  the  Xj_ 
direction . 

G 2 3 ^ G31,  G32  = Shear  moduli  in  the  2-3,  3-1,  and  1-2 

planes,  respectively. 

The  principal  material  coordinates  of  the  lamina  (x^, 
x2  , x^)  often  do  not  coincide  with  global  coordinates  (x,  y, 
z)  that  are  geometrically  natural  to  the  solution  of  the 
problem.  The  material  and  global  coordinates  are  related  by 


Figure  2.1.  Unidirectionally  reinforced  lamina 
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an  in-plane  rotation  as  shown  in  Figure  2.3,  where  0 is  the 
angle  measured  from  the  global  axis  x to  the  fiber  direction 
Xi . The  stress  and  strain  components  in  (x^,  x2 , x3  ) 
coordinates  are  related  to  the  stress  and  strain  components 
in  (x,  y,  z)  coordinates  by  the  tensor  transformation 
equations  [60].  The  transformations  are  commonly  written  as 


/ "s 

' s. 

al 

ax 

ct2 

ay 

a3 

°z 

= [ T ] 

0 4 

Tyz 

a5 

Txz 

JxY/ 

£1 

£ X 

e 2 

£y 

e 3 

Y23/2 

= [ T ] 

£ Z 

Yyz/2 

y31/2 

Yxz/2 

y12/2 

s / 

Y /2 

IJxy/  J 

and 


[ T 


o . ? 

cos  9 sin  9 0 

sin20  cos2g  0 
0 0 1 

0 0 0 

0 0 0 


-cosgsing  cosgsing  0 


0 0 2singcos9V 

0 0 -2sin0cosg 

0 0 0 

cosg  -sing  0 
sing  cosg  0 

0 0 cos^g-sin^g 


2.6 


23 


x3 


Figure  2.3.  Positive  rotation  of  principal  material 
axes  from  arbitrary  xy  axes 
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However,  if  the  matrix 
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due  to  Reuter  [61]  is  introduced,  then 
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By  combining  the  stress  and  strain  transformations  of 
Equations  2.4  and  2.5  along  with  Reuter's  matrix.  Equation 
, and  incorporating  the  strain-stress  relations  in 


2 . 7 
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principal  coordinates,  the  strain-stress  relations  in  the 
global  coordinates  can  be  established  in  the  form 
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where  superscript  -1  denotes  the  matrix  inverse.  By 
manipulating  the  matrix  multiplications  of  Equation  2.10, 
the  strain-stress  relations  for  an  arbitrarily  oriented 
lamina  with  respect  to  the  global  coordinates  (x,  y,  z) 

shown  in  Figure  2.3  are  obtained  as 
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Recall  that 

the  S ^ ^ are 

defined 

in  terms  of 

the 

engineering  constants  in  Equation  2.3.  The  transformed 
compliances  Sij  with  respect  to  the  global  coordinates  will 
be  used  in  the  development  of  the  hybrid  stress  finite 
element  method. 


2.2  Derivation  of  Stiffness  and  Mass  Matrices 
for  Hybrid  Stress  Models 

2.2.1  Introduction 

In  view  of  the  increasing  interest  in  using  composite 
materials  for  aerospace  structures,  the  analysis  of 
multilayered  plates  made  of  fiber-reinforced  composite 
materials  becomes  essential.  The  classical  lamination  theory 
( CLT ) , which  does  not  include  transverse  shear  effects,  can 
provide  reasonable  predictions  only  for  relatively  thin 
plates.  As  transverse  shear  effects  become  an  important  role 
in  laminated  plates  with  low  span-to-depth  ratios,  CLT  leads 
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to  a very  poor  description  of  laminate  response.  While 
analytical  solutions  [15,  16,  17]  are  usually  restricted  to 

problems  with  simple  loading  and  boundary  conditions,  the 
finite  element  method  is  found  attractive  in  dealing  with 
complicated  problems.  Laminated  plate  elements  based  on  the 
assumed  displacement  model  and  including  transverse  shear 
effects  have  been  developed  by  several  investigators  [5,  7, 

19,  21,  62].  It  has  been  noticed  that  there  always  exist 

some  limitations  and  difficulties  in  these  finite  element 

models.  The  above  mentioned  difficulties  can  be  overcome  by 
using  a hybrid  stress  finite  element  model.  Briefly,  the 

hybrid  stress  model  is  based  on  the  principle  of  minimum 
complementary  energy.  An  optimum  choice  of  the  number  of  the 
assumed  stress  modes  can  be  made,  which  gives  greater 
flexibility  in  the  descriptions  of  the  stress  field.  All  six 
stress  components  are  assumed  independently  within  each  ply 
element  and  satisfy  the  equations  of  equilibrium.  The 
formulation  of  the  hybrid  stress  model  for  both  a static 
case  and  a dynamic  case  will  be  presented  in  the  next 

section . 

2.2.2  Static  Case 

Based  on  the  development  of  Reissner  [63]  and  reference 
[22],  the  hybrid  functional  can  be  written  as 

"c  = l ^vn(aijeij  “ sijkl°ijakl/2)dv 
+ IavnTi(Qi-ui)ds  ' JsanfiGids] 


2.13 
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where 


a^j  = stress  tensor 
e^j  = strain  tensor 

= elastic  compliance  tensor 
= displacement  within  an  element 


u. 


= displacement  on  the  boundary  of  an  element 


vn  = volume  of  the  nth  element 

9vn  = boundary  of  the  nth  element 

San  = surface  of  the  nth  element  over  which 
traction  is  prescribed 

Tj_  = force  on  9vn 

= prescribed  boundary  traction 

The  Uj_  are  the  same  for  two  adjacent  elements  over  their 

common  boundaries  and  are  equal  to  prescribed  displacements 

u^  over  the  boundaries  where  the  displacements  are 

prescribed.  The  Reissner's  theorem  asserts  the  following: 

In  the  small-displacement  theory  of  elasticity,  the 
equilibrium  state  of  a body  is  such  that  6 tt c = 0 for 
arbitrary  variations  of  u^  and  a^j. 

The  Euler  equations  of  Equation  2.13  are 


(ui,  j + uj,i)//2 
°ij,j  = 0 
Ti  = °ijnj 


si jklakl 
in  v 


n 


2.14 


ui  = ui 


on  3v 


n 


where  nj  is  the  direction  cosine  of  a unit  vector  normal  to 
9vn  and 


T • = T • 

± l xi 


on  so 


n 


29 


In  general,  the  stresses  aj_j  and  the  displacements  Uj_ 
in  vR  and  the  displacements  u^  and  forces  Tj_  on  boundary  3vn 
may  be  independently  assumed.  It  will  lead  to  many  unknown 
variables,  which  may  or  may  not  be  useful  in  practice  from 
the  point  of  view  of  accuracy  and  efficiency.  To  reduce  the 
number  of  unknown  variables,  some  of  the  Euler  equations  are 
chosen  to  be  satisfied  identically. 

The  hybrid  functional  for  multilayered  elements  can  be 
expressed  in  matrix  notation  as 

TTC  = E [ J v (tTe  ■ aTSo/2)dv  - Jsa  T Tuds 
n n n 

+ J3V  TT(u  - u)ds  ] 2.15 

n 

where  superscript,  T,  denotes  the  matrix  transpose,  a is  the 
stress  column  matrix,  e is  the  strain  column  matrix  computed 
from  the  displacements  u,  S is  the  reduced  transformed 
compliance  matrix  as  described  in  Equations  2.11  and  2.12,  T 
is  the  column  matrix  of  surface  tractions,  T is  the  column 
matrix  of  prescribed  surface  tractions , and  u is  the  column 
matrix  of  boundary  displacements. 

The  stress  field  can  be  assumed  in  terms  of  unknown 
parameters  3 

a = P 3 2.16 

where  P is  a matrix  consisting  of  polynomial  terms  of 
coordinates  (x,  y,  z).  Displacements  u can  be  expressed  in 
terms  of  nodal  displacements  q 


u = N q 


2.17 
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where  N are  the  shape  functions  in  terms  of  the  natural 
coordinates  ( £,  n,  ?).  The  displacements  on  the  element 
boundary  can  be  described  as 

u = L q 2.18 

where  the  terms  in  matrix  L contain  coordinates  on  the 
surface.  The  strain  column  matrix  e can  be  calculated  from 
Equation  2.17  by  using  the  linear  strain-displacement 
relations 


e = Du  = DNq  = Bq  2.19 

Substituting  Equations  2.16-2.19  into  Equation  2.15  yields 

TTr  = I (BTGq  - BTHB/2  - qTQ ) 2.20 

° n 


where 

G = Jv  PTB  dv 

H = PTSP  dv  2.21 

J vn 

Q = I so  LtT  ds 
bUn 

For  simplicity,  the  above  formulations  are  written  on 
the  laminate  element  basis.  In  practice,  the  displacements 
and  stresses  are  assumed  within  a lamina  and  expressed  as 


oi  = p^i 
u^  = N^q^ 


2.22 


where  superscript,  i,  denotes  the  lamina  number  i.  It  is 
also  noted  that  the  stress  field  is  required  to  satisfy  the 
homogeneous  equilibrium  equations  identically.  Two 
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approaches  can  be  applied  in  the  development  of  a computer 
code  to  establish  the  relations  between  lamina  and  laminate. 
One  is  incorporating  the  interlaminar  surface  continuity 
conditions  (displacement  continuity,  interlaminar  traction 
continuity,  and  upper/lower  laminate  surface  traction-free 
condition)  to  relate  the  lamina  stress  parameters,  g1,  and 
nodal  displacements,  q1,  with  laminate  stress  parameters,  g, 
and  nodal  displacements,  q,  respectively,  and  then  using 
Equation  2.21  to  compute  matrices  G,  H,  and  Q.  The  other  is 
using  the  same  formulae  to  calculate  the  lamina  matrices, 
G1 , H1 , and  Q1  first,  then  with  proper  procedures  to 
assemble  the  lamina  matrices  into  laminate  matrices,  G,  H, 
and  Q.  The  details  of  stress  and  displacement  fields 
assumptions  and  interlaminar  surface  continuity  conditions 
will  be  described  in  the  next  section. 

Taking  variation  of  the  functional  ttc  (Equation  2.20) 
with  respect  to  g,  one  obtains 

G q - H g = 0 2.23 

and 

g = H- 1 G q 2.24 

By  substituting  g into  Equation  2.20,  the  hybrid  functional 
can  be  expressed  as 

tt  = Z ( qTKq/2  - qTQ ) 2.2  5 

° n 


where  Q is  the  static  loading  vector  (nodal  force  vector);  K 


32 


is  the  hybrid  element  stiffness  matrix,  defined  as 


K = GTH  1G 


2.26 


Taking  variation  of  Equation  2.25  with  respect  to  q yields 


Equation  2.27  is  the  typical  equation  for  the  finite  element 
method  in  static  problems.  Once  the  nodal  displacements  q 
have  been  computed,  stresses  and  strains  can  be  calculated 
by  using  Equations  2.24,  2.16,  and  2.11. 

2.2.3  Dynamic  Case 

A unified  approach  in  deriving  the  element  stiffness 
and  mass  matrices  for  dynamic  problems  based  on  the  hybrid 
stress  finite  element  model  has  been  presented  by  Tong  et 
al . [47].  The  functional  for  a modified  Reissner  principle 
[64,  65]  without  body  force  and  initial  stress  can  be 
defined  as 


where  u^  satisfy  the  prescribed  conditions  at  t = tj  and  t 
= t2  and  u^  = Uj_  over  the  boundaries  where  the  displacement 
is  prescribed  as  u^.  (*)  denotes  time  derivative,  (d/dt). 


independently  and  are  defined  the  same  as  in  static  case. 
Equation  2.13.  The  validity  of  Reissner 's  functional  is 


K q - Q = 0 


2.27 


2.28 


All  the  quantities,  a -j_  j , u^ , u^  , and  T^  , can  vary 
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based  on  the  existence  of  complementary  energy  and  the 
corresponding  constitutive  relations.  The  first  variation  of 
functional  is 


5tt 


ds 


R = l Vk  5ui(pV°ij,j)dv  + hvn  Sui(oijnj-Ti> 

+ Jv  - sljklakl1dv  + {u^.ds 


+ ST.(u.-u.)ds  - f 6u.T.ds]dt 

J 9 v 11  1 Jsa  ii 


2.29 


n n 

From  6ttr  = 0,  the  following  Euler  equations  can  be  obtained: 


PUi  - Oij;j  . 0 

(ui,j  + uj,i)/2  = sijkl  akl 


in  v 


n 


2.30 

on  3 vn 


The  vector  Tj_  is  continuous  at  the  common  boundary  of  any 
two  adjacent  elements  and 


Tj.  = T±  on  san 

If  the  a j_  j , Tj_ , u^,  and  u^  are  assumed  in  terms  of  unknown 
parameters  independently  in  the  finite  element  formulation, 
it  will  lead  to  many  unknown  variables.  To  reduce  the 
complexity,  a scheme  similar  to  the  hybrid  mixed  model 
suggested  by  Tong  et  al . [47]  will  be  employed.  The  first 
equation  (Equation  2.30)  is  not  actually  used,  since,  in  the 
finite  element  formulation,  all  the  mass  is  lumped  at  the 
nodes,  so  the  stress  in  the  elements  between  is  assumed  to 
satisfy  equilibrium  in  order  to  develop  the  stiffness 
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matrix.  Then  another  application  of  the  variational 
principle  will  be  used  to  develop  finite  element  equations 
of  motion  in  Equation  2.43. 

The  stresses 

separated  into  two  parts  as 

a . 


Oj_j  and  surface  tractions 


can  be 


i]  = °ij  + °ij 


m __  mH  i mF 

i Ai  ■Li 


2.31 


where 


13.1  >*  0 


T • H = a — n • 
1i  i]‘] 


m vn 
on  3 v 


2.32 


n 


Substituting  Equation  2.31  into  Equation  2.28  and  applying 
the  divergence  theorem  and  Equation  2.32,  transforms 
Equation  2.28  into 

= E f 2 [J  (a?.e..  - S.  ..  .o.  -a.  ,/2  -pu?/2)dv 
R n J t^  J vn  ij  ii  ljkl  ii  kl  M r 


+ J .e . -dv  + f 
Jvn  ii  ii  J 


- f T.u.ds]dt 
J sa  ii 
n 


„ T . (u. -u. ) ds 
3 v i i l 
n 


2.33 


In  practice,  one  may  make  the  following  assumption  to 
simplify  the  formulation  of  a hybrid  finite  element  model 
and  save  the  computation  time  [46]: 

ai?  = 0 2.34 

Substituting  Equation  2.34  into  Equation  2.33  yields 

”R  = 5 fq  [^vn<°ijeij  - Sijklaijakl//2  - P“2i/2)dv 

+ hv  Ti(ui-Ui)ds  - JSa  tAds!dt 


2.35 
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The  stress,  strain,  and  displacement  fields  are 
independently  assumed  as 

aH  = P 6 
Tf  = R a 
u = N q 

2.36 

u = N q 
u = L q 

e = Du  = DNq  = Bq 

The  assumptions  of  Equations  2.34  and  2.36  make  formulation 
approximate,  since  the  assumed  stress  field  does  not  satisfy 
the  equilibrium  equations  exactly;  instead  it  only  satisfies 
the  homogeneous  equilibrium  equations.  Substituting  Equation 
2.36  into  Equation  2.35,  transforms  the  hybrid  functional  ttr 
into  the  form 

u = Zjt2[BTGq  - BTHB/2  - qTMq/2  + aT  (G^q-G^)  - qTQ]dt  2.37 
n '"l 

where 

H = / v PTSP  dv 
n 

M = J v PNTN  dv 
n 


G = / v PTB  dv  2.38 


n 

" ^v 

RTL  ds 

n 

= ^9  v 

RTN  ds 

n 

- 

LTT  ds 

n 

functional  irR  with  respect  to  a and  3 yields 


The  variation  of 


G^q  G2q  — 0 
G q - H 3 = 0 


2.39 
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and 


G1  - g2 
3 = H_1G  q 


2.40 


A substitution  of  Equation  2.40  into  Equation  2.37  reduces 
the  functional  ttr  to  the  form 


TTp  = E /.  [q  Kq/2  - q Mq/2  - qT0]dt 
n n '~i 


2.41 


where 


K = GTH  1G 

M = { v pNTN  dv 
n 

Taking  the  first  variation  of  ttr  of  Equation  2.41 
respect  to  q yields 


2.42 


with 


K q + M q = Q 2.43 

Equation  2.43  are  the  equations  of  motion  for  the  finite 
element  method,  where  K,  M,  and  Q are  the  element  stiffness 
matrix,  mass  matrix,  and  the  load  vector,  respectively. 


2.3  Stress  and  Displacement  Interpolations 

for  a Lamina 

2.3.1  Introduction 

The  finite  element  analysis  is  based  on  the 

discretization  of  the  continuum  into  finite  subregions.  In 
the  conventional  displacement  method,  the  displacements  are 
chosen  as  the  prime  variables  and  the  stresses  are 

determined  from  the  calculated  displacement  field.  The 
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displacement  field  is  interpolated  in  terms  of  nodal 
displacements  through  some  interpolation  functions  also 
called  shape  functions.  The  choice  and  determination  of  the 
shape  functions  have  been  covered  by  most  of  the  finite 
element  method  texts  [66,  67,  68,  69,  70]. 

In  the  hybrid  stress  model,  both  stresses  and 
displacements  are  employed  as  variables  simultaneously.  The 
displacement  field  is  interpolated  through  shape  functions 
and  the  stress  field  is  interpolated  through  assumed  stress 
polynomials  with  unknown  stress  parameters.  Greater 
flexibility  is  achieved  in  the  stress  distributions  by  using 
the  hybrid  stress  model.  Punch  and  Atluri  [11]  have 
suggested  a systematic  procedure  to  develop  the  hybrid 
stress  model. 

Four-node  and  eight-node  hybrid  stress  laminated  plate 
elements  will  be  presented  in  this  section.  The  assumed 
stress  field  within  each  lamina  with  48  unknown  stress 
parameters  and  55  unknown  stress  parameters  will  also  be 
discussed . 

2.3.2  Displacement  Field  Interpolation 

Consider  the  multilayered  plate  element  composed  of  N 
perfectly  bonded  orthotropic  layers  in  which  the  various 
axes  of  elastic  symmetry  are  parallel  to  the  laminate  axes 
x,  y,  and  z as  shown  in  Figure  2.4,  where  z corresponds  to 
the  transverse  direction  of  the  laminate.  The  layers  are 
numbered  from  bottom  to  top  of  the  laminate.  The  interface 
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coordinates  are  denoted  as  z = h-^ , h2 , , hN+^_  . The 

thickness  and  orientation  of  composite  can  be  varied  from 
lamina  to  lamina.  To  achieve  the  general  applicability,  an 
isoparametric  finite  element  formulation  is  adopted.  The 
basic  idea  of  the  isoparametric  element  is  to  express  the 
element  coordinates  and  element  displacements  in  the  form  of 
interpolations  of  the  natural  coordinate  system  of  the 
element . 

For  the  present  multilayered  isoparametric  plate 
element,  the  coordinates  x and  y can  be  interpolated  as 


m 

X = .t  N -;X  • 

3=1  3 3 


y 


m 

j!lNjYj 


2.44 


where  x and  y are  the  in-plane  coordinates  at  any  point  of 
the  element,  and  x j , yj  (j=l,....,  m)  are  the  coordinates  of 
the  m nodal  points.  The  shape  functions  Nj  are  defined  in 
terms  of  the  natural  coordinates  ( £ , g ) of  the  element, 
£ and  n varying  from  -1  to  +1.  The  transverse  coordinate  z 
can  be  related  to  the  natural  coordinate  C in  the  form  as 


z = t(hi+hi+1)  + ?(hi+1-hi)]/2  2.45 

where  the  natural  coordinate  C within  each  lamina  also 
varies  from  -1  to  +1. 

The  displacements  u,  v,  and  w will  be  assumed  linearly 
through  the  thickness  of  each  lamina,  a manner  similar  to 
Mau  et  al.  [23]  and  Spilker  [32],  except  that  in  their 
assumptions  the  transverse  displacement  w was  constant 
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through  the  thickness  of  the  laminate.  In  a typical  lamina 

i,  the  displacement  field  is  related  to  the  natural 
coordinates  by 

m . . 

u1  =jI1Nj(^n.)[(l-c)u^/2  + (l  + ?)Uj1  + 1/2] 

vi  =jS1Nj(?,n)C(1"S)'Vi/2  + ( 1 + C ) v ji+1/2  ] 2.46 

w1  =j?  N j ( l ,n ) [ ( I"? ) Wj/2  + ( l+£ )Wj1+1/2 ] 

where  m is  the  number  of  element  nodes.  At  particular  node 

j,  the  degrees  of  freedom  are  shown  in  Figure  2.5.  Thus,  the 
total  number  of  degrees  of  freedom  per  node  is  3(N+1)  and 
the  total  number  of  degrees  of  freedom  per  element  is 
3m( N+l ) . 

For  the  four-node  isoparametric  plate  element  shown  in 
Figure  2.6,  the  shape  function  Nj  ( £ , ti ) corresponding  to  a 
specific  nodal  point  j is  expressed  as 

Nx  (£  , n ) = ( 1~£ ) ( l-n )/4 
Nt(^,  n ) = ( 1+5 )( i-n ) /4 

2.47 

n3(£,  n ) = ( 1+5  ) ( l+n.) /4 
n4(5,  n.)  = (1-5)  (l+n) /4 

In  the  eight-node  isoparametric  plate  element  shown  in 
Figure  2.7,  the  eight  shape  functions  are  described  as 

nx(£,  tv)  = (1-5)  (l-n)  (-l-5-n.)/4 
n2(5,  n)  = d~52) (l-n)/2 
n3(5,  n)  = (i+5)d-n)(-i+5-n)/4 
n4  (5  f n)  = d+5  ) d-n2)/2 


2.48 


40 


n 


(a)  Top  view 


Figure  2.4.  Geometry  and  layer  numbering  conventions  for 
the  isoparametric  multilayered  plate  element 


N+i 

— Uj 

N N N 
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3 3 3 

-V  V wj 
2 2 2 
— u . , v. , w. 

J J J 
l 1 ~ 1 

V-V  wj 

Nodal  degrees  of  freedom  for  the  displacement 
field 


Figure  2.5. 
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Figure  2.6.  A linear  (4-node)  isoparametric  element 


Figure  2.7.  A quadratic  (8-node)  isoparametric  element 


42 


N5(£,n)  = (1+5)  (l  + n)  (-l  + 5+n)/4 
n6(5  , n ) = ( 1-52)  ( i+n)/2 
n7(5,  n)  = (1- 5)  ( 1+n)  (-1- 5+n)/4 
n8(5  , n ) = (1-5)  ( 1 - n2 ) / 2 

Note  that  each  shape  function  has  the  value  of  one  at  its 
nodal  point  and  is  equal  to  zero  at  all  the  other  nodal 
points  in  both  cases. 

2.3.3  Stress  Field  Interpolation 

As  previously  mentioned,  all  six  stresses  are  included 
in  the  present  hybrid  stress  model.  The  stresses  within  each 
lamina  are  required  to  satisfy  the  homogeneous  equilibrium 
equations 

aij,j=0  2'49 

The  stress  distributions  within  each  lamina  are  interpolated 
in  terms  of  stress  polynomials  with  unknown  stress 
parameters,  8.  Various  choices  for  the  stress  interpolation 
polynomials  can  be  made  with  different  numbers  of  stress 
parameters,  depending  upon  the  structures  involved  and 
problem  requirements.  A necessary  but  not  sufficient 
condition  to  guarantee  a unique  solution  (i.e.  positive- 
definite  assembled  stiffness  matrix)  is  that  the  number  of 
stress  parameters  be  greater  than  or  equal  to  the  number  of 
generalized  nodal  displacements  minus  the  number  of  rigid 
body  modes  (equal  to  six  for  the  present  element)  [29].  Each 
extra  term  included  in  the  stress  assumption  will  add  more 
stiffness  [71]. 
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Two  versions  of  stress  field  within  each  lamina  will  be 
examined.  They  are  given  in  appendices  A and  B,  one  with  48 
stress  parameters  and  the  other  with  55  stress  parameters. 
The  in-plane  stresses  ax , Oy , and  rxy  are  assumed  to  vary 
linearly  in  the  transverse  direction  ? within  each  lamina. 
From  homogeneous  equilibrium  equations.  Equation  2.49,  t _ 
and  TyZ  will  be  of  order  £2  and  az  will  be  of  order 
within  the  lamina. 

For  the  present  multilayered  plate  element,  interface 
traction  continuity  and  laminate  upper/lower  surface 
traction  free  conditions  need  to  be  satisfied  exactly.  To 
ensure  these  conditions  being  satisfied  and  simplify  the 
writing  of  the  computer  code,  the  proper  numbering  of  the 
stress  parameters  is  required.  A scheme,  which  relates  the 
stress  parameters  in  each  lamina  to  the  laminate  element, 
suggested  by  Spilker  [32]  will  be  employed.  The  hybrid 
stress  model  with  48  stress  parameters  is  processed  as  an 
example.  The  application  to  the  55  stress  parameters  is 
similar . 

The  traction  free  conditions  on  the  upper  and  lower 
surfaces  of  the  laminate  in  the  present  multilayered  plate 
bending  element  are  described  as 

Txz  = Tvz  = 0 I toP  °f  the  laminate 

y 2.50 

az  = Txz  = Tyz  = 0 | bottom  of  the  laminate 

To  establish  the  interface  traction  continuity,  the 
following  condition  must  be  satisfied: 


top  of  the  ith  lamina 
bottom  of  the  (i+l)th  lamina 
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(c 


z ' 
z ' 


XT) 

xz ' yz ' 

T , T ) 

xz'  yz  ' 


From  the  stress  assumptions,  the  interlaminar  shear  and 
normal  stresses  at  upper  surface  (?  = 1)  of  lamina  i can  be 
expressed  in  the  form  as 


T xz  = (B37+338x+339Y+340xy)ti 
Tyz  = <341+342x+343Y+344xy)ti 
°z  = (B45+346x+347Y+B48xy)ti2 


and  at  the  lower  surface  (C=-l)  of  lamina  i yields 


Txz  = (^i+22x+33y+34xy)ti 

Tyz  = ( B5+B6x+37Y+38xy)ti  2'53 

az  = (B9+Biox+Blly+312XY)ti2 

where 

ti  = (hi+1-hi)/2  2.54 

To  satisfy  the  traction  free  conditions  at  lower  and  upper 
surface  of  the  laminate  requires  that 


^ = 0 j = 1,  2,  ,12 

= 0 j = 37,  38,  ,44 


2.55 


where  the  superscript  denotes  the  lamina  number  of  the 
laminate.  The  traction  continuity  at  the  interface  between 
lamina  i and  lamina  (i+1)  is  satisfied  by  demanding  that 
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B,i  + 1 = 


j 


j = 1,  2 


• • • • t 


8 


2.56 


j = 9,  10,  11,  12 


By  combining  the  laminate  surface  traction  free  and  the 
interface  traction  continuity  conditions.  Equations  2.55  and 
2.56,  with  the  stress  field  for  each  lamina,  some  stress 
parameters  can  be  eliminated.  Thus,  a set  of  stress 
parameters  of  the  multilayered  plate  element  is  formed.  In 
the  present  48-stress-parameter  model , the  total  number  of 
stress  parameters  for  the  N-layer  laminate  is  ( 48-12 )N-8. 
Similarly,  the  total  number  of  stress  parameters  will  be 
(55-14)N-10  for  the  55-stress-parameter  model. 

With  the  assumed  stress  and  displacement  fields , the 
element  stiffness  matrix  and  the  element  mass  matrix  can  be 
calculated.  The  static  and  dynamic  analysis  problems  can  be 
solved  by  using  the  proper  numerical  procedures  for  the 
solution  of  the  finite  element  equilibrium  equations.  The 
validity  of  the  hybrid  stress  finite  element  model  developed 
herein  can  be  determined  by  comparing  the  numerical  results 
with  the  existing  exact  analytical  solutions  [15,  16,  39]. 
Once  the  accuracy  of  the  finite  element  model  is 


established,  it  can  be  applied  to  general  laminate  analysis 
with  arbitrary  shape,  loading,  and  boundary  conditions,  for 
which  it  is  impossible  to  have  an  exact  analytical  solution. 


CHAPTER  3 


ANALYSIS  OF  STATIC  PROBLEMS 
3.1  Frontal  Solution  Method 

As  discussed  in  section  2.2.2,  for  a static  linear 
analysis  of  structures  and  solids,  the  finite  element 
formulation  yields  a set  of  linear  algebraic  equations  of 
the  form 


K q = Q 2.27 

where  K is  the  stiffness  matrix,  q is  the  nodal  displacement 
vector,  and  Q is  the  load  (nodal  force)  vector  of  the  finite 
element  system. 

The  method  employed  for  the  solution  of  the  equilibrium 
equations  is  a major  factor  influencing  the  efficiency  of 
any  finite  element  program.  The  techniques  which  are 
available  can  be  classified  as  either  direct  solution  method 
or  iterative  method.  For  the  present  study,  a specific 
solution  scheme  called  the  frontal  solution  method  will  be 
applied.  In  principle,  the  frontal  solution  method  is  a 
direct  elimination  process,  which  was  originated  by  Irons 
[72]  in  a 1970  paper  that  established  the  method  firmly  in 
finite  element  practice. 
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The  technique  can  be  considered  as  a particular 
technique  for  first  assembling  finite  element  stiffness 
matrices  and  nodal  force  vectors  into  a global  stiffness 
matrix  and  load  vector  and  then  solving  for  the  unknown 
displacements  by  means  of  a Gaussian  elimination  and 
backsubstitution  process.  It  is  designed  to  minimize  core 
storage  requirements,  the  number  of  arithmetic  operations 
and  the  use  of  peripheral  equipment.  The  disadvantage  is 
that  additional  out-of-core  operations  are  required,  which 
decrease  the  efficiency  of  the  method  by  a great  amount.  For 
the  details  of  the  frontal  solution  method,  a text  written 
by  Hinton  and  Owen  can  be  consulted  [73]. 

3.2  Example  Problems  and  Numerical  Results 
Two  types  of  example  problems  are  solved  by  the  hybrid 
stress  finite  element  method;  one  is  cylindrical  bending  of 
a simply  supported  long  strip  and  the  other  is  bending  of  a 
simply  supported  rectangular  laminated  composite  plate.  The 
accuracy  of  this  method  is  verified  by  comparing  the  results 
with  the  existing  three-dimensional  elasticity  exact 
analytical  solutions  by  Pagano  [15,  16]. 

3.2.1  Cylindrical  Bending  of  a Simply  Supported  Long  Strip 

The  laminated  plate  considered  in  this  example  consists 
of  layers  of  unidirectional  fibrous  composite  material.  The 
laminate  is  infinitely  long  in  the  y direction  and  simply 
supported  along  the  two  long  edges.  Sinusoidally  distributed 
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transverse  loading  q0sin(^x/L)  is  applied  on  the  top  surface 
of  the  laminate  as  shown  in  Figure  3.1.  Three  different 
laminate  configurations  are  considered, 

(1)  a two-layer  cross-ply  laminate  (0/90)  with  layers 
of  equal  thickness, 

(2)  a three-layer  symmetric  cross-ply  laminate  (0/90/0) 
with  layers  of  equal  thickness, 

(3)  a four-layer  antisymmetric  cross-ply  laminate  (90 
/0/90/0)  with  layers  of  equal  thickness. 

In  all  the  problems,  the  material  properties  for  each 
lamina  are 


'L  = 

174.6 

GPa 

( 

25xl06 

psi 

'T  = 

7 GPa 

(10 

6 

psi  ) 

;LT 

= 3.5 

GPa 

(0 

. 5xl06 

psi 

'TT 

= 1.4 

GPa 

(0 

. 2xl06 

ps  i 

LT 

= VTT 

= 0. 

25 

3.1 


where  L denotes  the  direction  parallel  to 
the  transverse  direction.  For  each  case, 
thickness  ratios  (S)  are  investigated.  The 
will  be  presented  in  terms  of  normalized 
defined  as 


the  fibers  and  T 
various  span-to- 
numerical  results 
values  which  are 


ax  = ax(L/2,z)/q0 
cz  = cz(L/2,z)/q0 

Txz  “ Txz^'z^  /^o 


u = ETu(o,z)/hqQ 


3.2 


49 


(b)  Side  view 


Figure  3.1.  Cylindrical  bending  of  a simply  supported 
long  strip  under  sinusoidal  loading 
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w = 100ETh^w( L/2 , z ) /qQL4 
S = L/h  z = z/h 

The  normalized  central  transverse  deflections  w with 
respect  to  S for  the  three  cases  of  cylindrical  bending  are 
shown  in  Table  3.1,  where  the  surface  number  indicates  the 
location  of  each  surface  of  laminate  starting  from  bottom  to 
top  of  the  laminate  (see  Figure  3.1).  The  elasticity 
solutions  are  evaluated  by  using  the  same  procedures  as 
given  in  Pagano  [15].  The  classical  lamination  theory 
results,  which  are  independent  on  the  span-to-thickness 
ratio,  are  also  presented  in  case  2.  The  results  predicted 
by  the  hybrid  stress  model  are  in  good  agreement  with  the 
exact  elasticity  solutions  while  the  CLT  underestimates  the 
transverse  deflection.  The  results  evaluated  by  CLT  are 
accurate  only  for  the  thin  plate  (i.e.  S is  larger  than  50). 
From  Table  3.1,  it  is  noted  that  the  assumption  of  constant 
transverse  deflection  through  the  laminate  thickness  can  not 
exactly  describe  the  response  for  the  thick  laminate  (for 
example,  S equal  to  4).  The  difference  of  transverse 
deflection  between  the  top  surface  and  bottom  surface  of  the 
laminate  is  varied  from  3%  to  7%  depending  on  the  geometry 
configuration  of  the  laminate.  For  case  2,  the  normalized 
horizontal  displacement  u and  normalized  stress 
distributions  qx,  t xz > and  gz  through  the  laminate  thickness 
for  S = 4 and  S = 10  are  shown  in  Figures  3.2  and  3.3, 
respectively.  Excellent  agreements  are  found  for  the  hybrid 
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Table  3.1.  Effects  of  span-to-thickness  ratio  S on  the 

normalized  central  transverse  deflections  for 
the  cylindrical  bending  problem 


Case 

no. 

S = 

★ 

Pagano 

4 

Present 

S = 
Pagano 

10 

Present 

S = 
Pagano 

50 

Present 

CLT 

3 

4 

.760 

4.773 

2.9  42 

2.9  44 

2.634 

2.632 

1 

2 

4 

.695 

4.703 

2.954 

2.954 

2.634 

2.633 

1 

4 

.623 

4.595 

2.951 

2.9  52 

2.634 

2.633 

4 

3 

.023 

3.022 

0.9  34 

0.9  33 

0.527 

0.527 

0, 

.5096 

9 

3 

2 

.925 

2.931 

0.933 

0.932 

0.527 

0.527 

0, 

.5096 

2 

2 

.864 

2.868 

0.931 

0.931 

0.527 

0.527 

0. 

.5096 

1 

2 

.839 

2.849 

0.929 

0.927 

0.527 

0.527 

0. 

.5096 

5 

4 

.299 

4.295 

1.661 

1.661 

1.146 

1.147 

4 

4 

.229 

4.230 

1.661 

1.663 

1.146 

1.147 

3 

3 

4 

.181 

4.184 

1.660 

1.663 

1.146 

1.147 

2 

4 

.156 

4.162 

1.658 

1.661 

1.146 

1.147 

1 

4 

.115 

4.117 

1.654 

1.656 

1.146 

1.147 

Pagano:  Elasticity  exact  solution  [15] 

Present:  Hybrid  stress  finite  element  method 

CLT:  Classical  lamination  theory 


★ 


no.  : 


Surface  number 
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(a)  Normal  displacement 


z 


(b)  Normal  stress 

Figure  3.2.  Through-thickness  normalized  displacement  and 
stress  distributions  for  3=4,  case  2 
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z 


(c)  Transverse  shear  stress 


(d)  Transverse  normal  stress 


Figure  3.2.  (Cont.) 
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(a)  Normal  displacement 


(b)  Normal  stress 

Figure  3.3.  Through-thickness  normalized  displacement  and 
stress  distributions  for  S=10,  case  2 
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z 


(c)  Transverse  shear  stress 


Interface 


a 


z 


(d)  Transverse  normal  stress 


Figure  3.3.  (Cont.) 
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stress  model  and  the  exact  elasticity  solution.  The 
deviations  between  CLT  and  elasticity  solution  are  very 
significant  for  the  case  S equals  4.  The  CLT  results 
converge  to  exact  elasticity  solutions  with  increasing  span- 
to-thickness  ratio.  Figures  3.2a  and  3.3a  indicate  that  the 
assumption  of  linear  displacement  variation  within  each 
lamina  is  reasonable.  In  Figure  3.4,  the  normalized 
horizontal  displacement  and  normalized  stresses  for  S = 4 of 
case  3 are  plotted.  The  accuracy  of  the  hybrid  stress  model 
with  increasing  layer  number  is  confirmed  again.  Note  that 
the  transverse  normal  stress  which  cannot  be  predicted  by 
[23]  is  also  presented  in  Figures  3. 2d,  3.3d,  and  3.4d. 

3.2.2  Bending  of  a Simply  Supported 
Rectangular  Laminated  Plate 

This  type  of  problem  has  been  investigated  by  Pagano 
[16];  the  configuration  of  the  laminated  plate  is  shown  in 
Figure  3.5.  To  validate  the  accuracy  of  the  present  method, 
the  results  are  compared  with  Pagano  [16].  Two  different 
geometries  of  the  laminated  plate  with  various  S ( = a/h ) 
will  be  under  investigation,  i.e., 

(1)  a symmetric  3-layer  cross-ply  (0/90/0)  laminated 
plate  with  a = b and  layers  of  equal  thickness, 

(2)  the  same  lamination  geometry  as  in  (1),  but  the 
laminate  is  of  rectangular  shape,  b = 3a . 

The  plates  are  subject  to  sinusoidally  distributed  loading. 
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(a)  Normal  displacement 


(b)  Normal  stress 


Figure  3.4. 


Through- thickness  normalized  displacement  and 
stress  distributions  for  S=4,  case  3 
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(c)  Transverse  shear  stress 


Interface 


a 


z 


(d)  Transverse  normal  stress 


Figure  3.4.  (Cont.) 
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(a)  Top  view 


(b)  Side  view 

Figure  3.5.  Bending  of  a simply  supported  rectangular 
laminated  plate  under  sinusoidal  loading 
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q(x,y)  = q0sin(TTx/a)sin(TTy/b)  3.3 

The  lamina  material  properties  are  the  same  as  in  the 
previous  section.  Equation  3.1.  The  predicted  results  are 
presented  in  Tables  3.2  and  3.3,  where  the  elasticity 
solution  [16],  higher  order  shear  deformation  theory 
solution  investigated  by  Reddy  [21],  and  CLT  solution  are 
also  included  for  comparison.  For  the  case  S = 4,  the 

maximum  tX2  does  not  occur  at  z = 0;  therefore  the  second 
entry  represents  the  maximum  value  with  the  coordinates  z 
being  given  in  parentheses.  The  normalized  quantities  are 
expressed  as 

(V  V fxy}  = 'ax'  V Txy)/c3os2 

^xz'  ^yz  ^ - ^Txz  ' Tyz^//c3o^ 

w = 100ETw/qohS4 

S = a/h 

From  the  results  summarized  in  Tables  3.2  and  3.3,  the 
excellent  accuracy  is  observed  for  the  hybrid  stress  finite 
element  method.  The  CLT  solution  is  accurate  only  for  thin 
plates;  the  disagreement  between  elasticity  solution  and  CLT 
solution  is  significant  when  the  span-to-thickness  ratio  S 
is  lower. 
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Table  3.2.  Normalized  deflection  and  stresses  in  a 3-layer 
(0/90/0)  square  laminate  (a  = b)  subject  to 
sinusoidal  loading 


- .a 
ax(2\ 

- .a 

ay  2, 

X T 

xz  vz 

- ,a 

x w (-=• 

xy  2 , 

s 

Source 

k +*L) 
2,  -2 

k +h) 
2,  -6J 

o 

Mtr 

o 

O 

o 

(0,0,  +|)  | 0 

.801 

.534 

.256 

-.0511 

Pagano 

.217 

-.755 

-.556 

. 282 ( . 27h ) 

.0505 

. 717 

.517 

.263 

-.0476 

4 

Present 

.221 

2.020 

-.679 

-.541 

. 284 ( . 27h ) 

.0485 

Reddy 

.7345 

.1832 

1.9218 

.285 

Pagano 

±.590 

-.288 

.357  .1228 

+.0289 

.285 

-.0283 

10 

Present 

±.580 

.367  .127 

.7546 

-.289 

.0284 

Reddy 

.5684 

.1033 

.7125 

Pagano 

±.552 

±.210 

.385  .0938 

+.0234 

20 

Present 

±.553 

±.210 

.395  .0998 

+.0233  .5170 

Reddy 

CLT 

±.539 

±.180 

.395  .0823 

+.0213 

Pagano : 
Present: 
Reddy: 
CLT : 


Elasticity  exact  solution  [16] 

Hybrid  stress  finite  element  method 
High  order  shear  deformation  theory  [21] 
Classical  lamination  theory 
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Table  3.3.  Normalized  deflection  and  stresses  in  a 3-layer 
(0/90/0)  rectangular  laminate  (b  = 3a)  subject 
to  sinusoidal  loading 


- /& 
° x 2 , 

- .a 

°y  2, 

T 

xz 

Tvz 

T 

xy 

»<f. 

s 

Source 

k +h) 
2,-2 

k +h) 
2,-6; 

(0,|,0)  (• 

t'0'0) 

(0,0, +|) 

!,o) 

1.14 

.109 

.351 

-.0269 

" 

Pagano 

-1.10 

-.119 

. 387 ( . 2 7h ) 

.0334 

.0281 

2.82 

1.077 

.108 

.360 

-.0254 

4 

Present 

-.975 

-.118 

. 390 ( . 27h ) 

.0326 

.0267 

2.828 

Reddy 

1.0356 

.1028 

.2724 

.0348 

.0263 

2.6411 

.726 

.0418 

-.0120 

Pagano 

-.725 

-.0435 

.420 

.0152 

.0123 

.919 

.709 

.0429 

-.0118 

10 

Present 

-.707 

-.0448 

.428 

.0151 

.0121 

.921 

Reddy 

.6924 

.0398 

.2859 

.0170 

.0115 

.8622 

.0294 

Pagano 

+ .650 

-.0299 

.434 

.0119 

+.0093 

.610 

.653 

.0298 

-.0092 

20 

Present 

-.646 

-.0304 

.450 

.0118 

.0093 

.611 

Reddy 

.6407 

.0289 

.2880 

.0139 

.0091 

.5937 

CLT 

+ .623 

+.0252 

.440 

.0108 

+.0083 

.503 

Pagano: 

Elasticity  exact  solution 

[16] 

Present:  Hybrid  stress  finite  element  method 


Reddy:  High  order  shear  deformation  theory  [21] 

CLT:  Classical  lamination  theory 
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3.3  Characteristics  of  Present 
Hybrid  Stress  Finite  Element  Model 

All  the  problems  in  previous  section  were  solved  by  the 
eight-node  isoparametric  element  with  55  stress  parameters 
for  each  individual  lamina.  Solutions  by  using  the  four-node 
isoparametric  element  with  48  stress  parameters  were  also 
obtained.  The  results  predicted  by  the  eight-node 
isoparametric  element  were  only  slightly  more  accurate  than 
those  predicted  by  the  four-node  isoparametric  element.  The 
maximum  difference  between  four-node  and  eight-node 
isoparametric  element  results  is  less  than  2 per  cent. 

In  the  cylindrical  bending  of  a long  strip  problem,  4 
elements,  8 elements,  12  elements,  and  16  elements  in  half 
span  have  been  adopted  to  compare  the  convergence  of  the 
results.  It  appears  that  the  convergence  rate  of  the 
displacement  field  is  faster  than  that  of  the  stress  field. 
For  all  the  span-to-thickness  ratios,  the  displacement  field 
can  be  predicted  accurately  by  only  4 elements  in  the  half 
span.  However,  in  the  stress  field,  12  elements  in  the  half 
span  were  required  for  the  thin  plate  case  (S  = 50)  while  4 
elements  were  enough  to  get  an  accurate  stress  field  for  the 
thick  plate  case  (S  = 4,  or  S = 10).  Therefore,  in  order  to 
obtain  an  accurate  stress  field  in  the  thin  plate,  a more 
refined  mesh  is  required.  The  same  tendency  was  also 
observed  in  the  bending  of  a simply  supported  rectangular 
laminated  plate. 

From  the  numerical  results  of  the  previous  section,  the 
accuracy  of  the  present  hybrid  stress  finite  element  method 
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has  been  established.  The  present  approach  is  sufficiently 
general  to  describe  the  laminate  elastic  response  and  can  be 
implemented  easily  to  attack  complicated  problems  which 
cannot  be  solved  exactly. 


CHAPTER  4 


DYNAMIC  ANALYSIS  OF  LAMINATED  PLATES 
4.1  Introduction 

The  derivation  of  element  stiffness  and  mass  matrices 
by  using  the  hybrid  stress  model  has  been  discussed  in 
section  2.2.3.  In  this  chapter  the  analyses  of  the  free 
vibration  and  transient  response  of  a laminated  plate  are 
presented.  Simply  supported  square  laminated  anisotropic 
plates  having  both  symmetric  and  anti-symmetric  laminations 
are  considered  for  the  natural  frequencies  analysis.  The 
effect  of  the  degree  of  orthotropy  of  the  individual  layer, 
the  effect  of  number  of  layers,  and  the  effect  of  length  to 
thickness  ratio  on  the  fundamental  frequency  are  examined. 
In  transient  analysis  , the  responses  of  simply  supported 
laminated  plates  subject  to  sinusoidal  loading  and  uniform 
loading  are  investigated.  For  the  solution  of  the  finite 
element  equations  of  motion  of  free  vibration  and  transient 
response  problems,  two  effective  solution  techniques,  the 
subspace  iteration  method  [68]  and  the  Newmark  direct 
integration  method  [74]  will  be  presented,  respectively. 
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4.2  Subspace  Iteration  Method 
The  free  vibration  finite  element  equations  of  motion 
with  damping  neglected  are 

M q + K q = 0 4.1 


where  K is  the  structure  stiffness  matrix  and  M is  the 
structure  mass  matrix.  Equation  4.1  can  be  solved  by 
expressing  the  field  variable  as 


q = 


4.2 


where  <f>  is  a nodal  vector  of  order  n,  t the  time  variable, 
and  a)  the  natural  frequency  of  vibration  of  the  plate  in  the 
mode  described  by  the  vector  <f> . Substituting  Equation  4.2 
into  Equation  4.1  yields  the  generalized  eigenvalue  problem 

K - «2  M <}>  = 0 4.3 


from  which  <P  and  m can  be  determined.  For  matrices  of 

dimension  n x n,  there  will  be  n eigensolut  ions  (m^2,  4>  ^ ) r 
2 2 

( oj2  , <j>2 ) ( “n  ' ^n ) • An  important  property  of  the 

eigenvectors  is  that  they  satisfy  the  orthogonality 

conditions,  i.e. 


bj_T  M 4>j 
<f>iT  K (frj 


and 


0 < co-^2  < 


< 


4 . 4 


There 


are  many  different  techniques  existing  for  the 
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solution  of  eigenvalue  problems.  Since  the  procedures  for 
the  eigenvalue  problems  are  time  consuming,  the  choice  of  an 
appropriate  and  effective  method  is  an  important  factor  for 
the  general  application,  especially  in  the  large  eigenvalue 
problem.  The  subspace  iteration  method  suggested  by  Bathe 
[68]  will  be  adopted  to  conduct  the  investigation.  This 
method  has  been  used  extensively  in  a number  of  general- 
purpose  finite  element  analysis  programs  and  has  proven 
cost-effective  and  reliable.  In  structure  analysis,  the 
lowest  few  eigenvalues  (natural  frequencies)  are  the  main 
concern  of  investigators.  The  basic  objective  in  the 
subspace  iteration  method  is  to  solve  for  the  p smallest 
eigenvalues  and  corresponding  eigenvectors,  which  satisfy 

K $ = M $ A 4.5 

where  $ = [ cj)1 , <j >2 , <|>  ] 

and  A is  a diagonal  matrix  of  and  the  eigenvector 

also  satisfies  the  orthogonality  conditions  (Equation  4.4). 

The  subspace  iteration  method  consists  of  three  steps 

[68  ] : 

1.  Establish  q starting  iteration  vectors;  q > p,  q = 
min(2p,  p+8)  is  a proper  selection,  where  p is  the 
number  of  eigenvalues  and  eigenvectors  to  be 
calculated . 

2.  Use  simultaneous  inverse  iteration  on  the  q vectors 
and  Ritz  analysis  to  extract  the  "best"  eigenvalue 
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and  eigenvector  approximations  from  the  q iteration 
vectors . 

For  k = 1,  2 , 

K Xk+1  = M Xk  4.6 

where  X^  is  the  starting  iteration  vector. 

Find  the  projections  of  the  operators  K and  M, 

Kk+i  - 5k+iT  K xk+1 

— m — 4.7 

\+l  = *k+l  M Xk+1 

Solve  for  the  eigensystem  of  the  projected  operators, 

Kk+1  Qk+1  = Mk+1  Qk+1  Ak+1  4-8 

Find  an  improved  approximation  to  the  eigenvectors, 

Xk+1  = xk+l  Qk+1  4.9 

As  k + =o 


Ak+1  - A and  xk+l  ^ $ 

3.  After  iteration  convergence,  use  the  Sturm  sequence 
check  to  verify  that  the  required  eigenvalues  and 
corresponding  eigenvectors  have  been  calculated. 
Reference  [68]  has  presented  very  detailed  descriptions 
about  the  subspace  iteration  method  and  is  a good  reference 
to  use  to  become  familiar  with  this  method. 
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4.3  Newmark  Direct  Integration  Method 
As  derived  in  section  2.2.3,  the  governing  finite 
element  equations  of  motion  for  a linear  dynamic  analysis 
are 

M q + K q = Q 2.43 

where  M and  K are  the  mass  and  stiffness  matrices,  q and  q 
are  the  nodal  displacement  and  acceleration  vector,  and  Q is 
the  load  (nodal  force)  vector  of  the  finite  element  system. 
In  finite  element  analysis,  there  exist  many  effective 
numerical  procedures  to  solve  the  linear  differential 
equations,  Equation  2.43.  Basically,  they  can  be  divided 
into  two  methods  of  solution:  the  direct  integration  method 

and  the  mode  superposition  method.  In  the  present  study,  the 
Newmark  direct  integration  method  [74]  will  be  followed  to 
integrate  Equation  2.43  step  by  step.  The  following 
assumptions  are  made  in  the  numerical  analysis: 

^n+1  = <3n  + [ ( l-S)qn+6qn+1]At  4.10 

Sn+l  = % + <3nAt  + 2 ~a  )%+a%+l ] At2  4.11 

where  At  is  the  time  step  size,  n is  the  step  number,  and 
the  parameters  6 and  a control  integration  accuracy  and 
stability.  At  time  tn+^  = (n+l)At,  the  finite  element 
equations  of  motion  (Equation  2.43)  are  described  as 

M %+i  + K qn+i  = Qn+i  4.i2 

Solving  from  Equation  4.11  for  qn+1  in  terms  of  q +1 , and 
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then  substituting  into  Equation  4.12  and  rearranging  the 
terms  transforms  the  equations  to  the  form 


where 


and 


k gn+i  = Qn+i 


4.13 


Qn+1  = Qn+1  + M(ao%+alc3n+a2'<3n) 


aQ  = 1/a  (At ) 2 

ax  = 1/a (At ) 4.15 

a2  = (1  - 2a) /2a 


Once  the  displacements  qn+1  at  time  step  n+1  are  known  from 
Equation  4.13,  the  velocities  and  accelerations  can  be 
computed  by  using  Equation  4.10  and  4.11  and  expressed  as 

%+l  = ao((3n+l  " " al%  " a2%  4‘16 

<5n+l  = <3n  + a3%  + a4%+l  4 • 17 

where 


a 3 = (1  - 5) (At) 
a4  = 6 ( At ) 


4 .18 


A special  scheme  originated  by  Newmark  with  6=0.5  and  a = 
0.25  is  used  here  to  integrate  the  equations  step  by  step. 
These  values  correspond  to  the  constant-average-acceleration 
method,  which  gives  an  unconditionally  stable  numerical 
scheme  [ 74 ] . 
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4.4  Numerical  Examples 

4.4.1  Free  Vibration  of  a Simply  Supported 
Laminated  Square  Plate 

Fundamental  frequencies  of  simply  supported  square 
orthotropic  laminates  having  both  symmetric  and  anti- 
symmetric laminations  with  respect  to  the  mid-plane  are 
considered.  Because  of  the  existing  of  biaxial  symmetry  in 
the  cross-ply  laminate,  only  a quadrant  of  the  laminate  is 
modeled  in  the  finite  element  analysis.  In  both  laminations, 
the  total  thickness  of  0°  laminae  is  equal  to  the  total 
thickness  of  the  90°  laminae  in  each  laminate.  The  boundary 
conditions  for  a quadrant  and  the  geometry  configurations 
for  both  laminations  are  shown  in  Figure  4.1.  The  material 
properties  of  each  lamina  are 

El/Et  = open 
Get/^t  = 0.6 

Gtt/Et  = 0.5  4.19 

Vrpip  0.25 
p ( density)  = 1.0 

The  effects  of  degree  of  orthotropy,  EL/ET,  of  the 
individual  layer  and  the  effect  of  number  of  layers  on  the 
fundamental  frequency  of  symmetric  and  anti-symmetric  cross- 
ply  laminates  have  been  examined.  Table  4.1  contains  the 
fundamental  frequencies  predicted  by  the  hybrid  stress 
finite  element  method  for  square  cross-ply  laminates  (a/h  = 
5).  The  results  are  compared  with  Noor 's  [39]  solutions  of 


u=0 


(a)  Simply  supported  boundary  conditions  for  a quadrant 
of  the  laminate 


z,w 


T 

h 

1 


a/2 


90 


_21L 


90_ 


x,u 


(b)  Anti-symmetric  lamination  of  the  laminate 


(c)  Symmetric  lamination  of  the  laminate 


Figure  4.1.  Boundary  conditions  and  the  geometry 
configurations  of  the  laminate 
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Table  4.1.  Effect  of  degree  of  orthotropy  of  the  individual 
layer  on  the  fundamental  frequency  of  simply 
supported  square  multilayered  composite  plates 
with  a/h  = 5 


No . of 

Source  Layers  3 

10 

el/et 

20 

30 

40 

Noor 

.25031 

.27938 

.30698 

.32705 

.34250 

Present 

.25148 

.28030 

.30768 

.32763 

.34301 

Reddy 

2 

.27955 

.31284 

.34020 

.36348 

CLT 

.30968 

.35422 

.39335 

.42884 

Noor 

.26182 

.32578 

.37622 

.40660 

.42719 

Present 

.26219 

.32621 

.37675 

.40719 

.42780 

Reddy 

4 

.32782 

.38506 

.42139 

.44686 

CLT 

.38877 

.49907 

.58900 

.66690 

Noor 

.26440 

.33657 

.39359 

.42783 

.45091 

Present 

.26458 

.33666 

.39359 

.42775 

.45077 

Reddy 

6 

.33621 

.39672 

.43419 

.46005 

CLT 

.40215 

.52234 

.61963 

.70359 

Noor 

Present 

Reddy 

CLT 

3 

.26474 

.26524 

.32841 

.33364 

.33087 

.41264 

.38241 

.38289 

.38105 

.54043 

.41089 

.41142 

.41087 

.64336 

.43006 

.43062 

.43149 

.73196 

Noor 

.26587 

.34089 

.39792 

.43140 

.45374 

Present 

.26608 

.34103 

.39803 

.43149 

.45380 

Reddy 

5 

.33988 

.39935 

.43502 

.45917 

CLT 

.41264 

.54043 

.64336 

.73196 

0)  =w  (ph2/ET)1//2 

Noor:  Elasticity  solution  [39] 

Present:  Hybrid  stress  finite  element  method 

Reddy:  Refined  shear  deformation  theory  [45] 

Classical  lamination  theory  [45] 


CLT : 
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the  three-dimensional  elasticity  theory  by  using  higher 
order  finite  difference  schemes.  Very  good  agreements  are 
observed  between  the  hybrid  stress  model  and  3-D  elasticity 
theory.  The  results  obtained  by  using  classical  lamination 
theory  and  refined  shear  deformation  theory  [45]  which  takes 
into  account  the  shear  deformation  and  rotary  inertia  are 
also  included  for  comparison.  The  CLT  overestimates  the 
fundamental  frequencies,  especially  when  the  degree  of 
orthotropy  is  greater.  Figure  4.2  shows  the  influence  of  the 
degree  of  orthotropy  of  the  individual  layer  on  the  accuracy 
of  the  fundamental  frequencies  obtained  by  the  CLT  for 
laminates  with  both  symmetric  and  anti-symmetric  lamination. 

The  effect  of  length-to-thicknes s ratio  on  the 
fundamental  frequency  for  different  EL/ET  of  the  two-layer 
anti-symmetric  and  three-layer  symmetric  cross-ply  laminates 
is  shown  in  Figure  4.3,  where  the  normalized  fundamental 
frequency  is  expressed  as 

5 = (o)2pa4/ETh2  ) X/2  4.20 

Severe  differences  in  the  normalized  fundamental  frequencies 
are  found  at  lower  length-to-thickness  ratios  (for  example 
a/h  < 20).  The  difference  is  dependent  on  the  degree  of 
orthotropy;  the  greater  value  of  EL/ET  produces  a larger 
difference.  When  a/h  is  higher  than  20,  there  is  no 
significant  change  in  the  normalized  fundamental  frequency 
with  increasing  length-to-thickness  ratio.  It  is  seen  that 
the  normalized  fundamental  frequency  will  approach  the  CLT 
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«c  = CLT  fundamental  frequency 

= Hybrid  stress  model  fundamental  frequency 

Figure  4.2.  Effect  of  Et/Et  on  the  accuracy  of  the 
fundamental  frequencies  obtained  by  CLT 
for  laminated  orthotropic  square  plates 
with  a/h=5 
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(a)  Two-layer  anti-symmetric  laminate 


(b)  Three-layer  symmetric  laminate 


Figure  4.3.  Effect  of  a/h  on  the  fundamental  frequencies 
of  square  cross-ply  laminates 
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result  asymptotically  for  higher  values  of  a/h  (thin 
laminated  plates). 

A hybrid  stress  model  with  a four-node  isoparametric 
element  and  55  assumed  stress  parameters  for  each  individual 
layer  was  used  in  the  present  investigation.  Considering  the 
convergence  rate  of  the  present  model , the  same  tendency  is 
observed  as  in  the  static  problem.  For  a thin  laminate,  a 
more  refined  mesh  is  required  to  obtain  an  accurate  result. 

In  conclusion,  the  results  of  the  present  study  show 
that  the  hybrid  stress  finite  element  model  can  predict  the 
free  vibration  behavior  of  an  orthotropic  laminate  very 
accurately;  the  validity  of  the  hybrid  stress  model  in 
dynamic  analysis  of  laminated  plates  is  confirmed.  Again, 
the  CLT  is  not  adequate  for  the  analysis  of  free  vibration 
of  thick  laminate  (a/h  < 20). 

4.4.2  Transient  Response  of  a Simply  Supported 
Laminated  Square  Plate 

The  transient  response  of  simply  supported  laminated 
plates  is  presented  in  this  section.  The  laminates  are 
subjected  to  suddenly  applied  sinusoidally  distributed  pulse 
loading , 

q ( x , y,  t)  = (q0sinTTx/a  s iniry/b ) H ( t ) 4.21 

where  H(t)  is  the  Heavyside  step  function.  The  following  two 
laminated  plates  are  considered: 

1.  A two-layer  anti-symmetric  cross-ply  (90/0)  square 
laminate  with  layers  of  equal  thickness. 
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2.  A three-layer  symmetric  cross-ply  (0/90/0)  square 
laminate  with  layers  of  equal  thickness. 

In  both  problems , the  same  material  properties  as  in  Putcha 
and  Reddy  [45]  are  employed  for  each  individual  layer, 

El  = 525  GPa 
Eip  = 21  GP  a 

Gj^ip  = Grp rp  = 10.5  GPa  4.22 

Vp^rp  — Vrprp  = 0 • 25 
p = 0 . 8 g/cm2 

Owing  to  the  biaxial  symmetry  of  the  laminate  geometry,  only 
one  quadrant  of  the  laminate  is  analyzed.  The  geometry 
configurations  and  boundary  conditions  of  the  finite  element 
model  are  shown  in  Figure  4.4. 

The  normalized  deflection  and  stresses  are  described  as 

w = 1000ETh^w/qoa4 

(ax,  a y)  = 10(ax,  ay)/qQS2 

4.23 

^ az ' Txz  ' Tyz  ^ - Txz'  Tyz^//<3o^ 

S = a/h  z = z/h 

In  the  present  study  qQ  is  taken  to  be  100  and  the  time  step 
size  is  equal  to  5 microseconds.  The  normalized  central 
transverse  deflection,  w(a/2,  a/2,  0),  as  a function  of  time 
for  a two-layer  anti-symmetric  simply  supported  cross-ply 
square  laminate  under  sinusoidal  loading  is  shown  in  Figure 
4.5.  Throughout  Figures  4.6,  4.7,  4.8,  and  4.9,  the  stresses 
with  respect  to  time  for  a two-layer  laminate  are  plotted. 
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y.v 


(a)  Top  view  and  boundary  conditions  of  laminate 


z 


T 

0 

1 h/2 

n 

90 

_h/2 

(b)  Two-layer  anti-symmetric  laminate 


z 


(c)  Three-layer  symmetric  laminate 


Figure  4.4.  Laminate  geometry  configurations  and  boundary 
conditions 
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Figure  4.5.  Center  deflection  versus  time  for  a 2-layer  cross-ply  (90/0)  simply 
supported  square  laminate  under  suddenly  applied  sinusoidal  loading 
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Figure  4.6.  Normal  stress  versus  time  for  a 2-layer  (90/0)  laminate  under  suddenly 
applied  sinusoidal  loading 
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Figure  4.7.  Normal  stress  versus  time  for  a 2-layer  (90/0)  laminate  under  suddenly 
applied  sinusoidal  loading 
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Figure  4.8.  Transverse  shear  stresses  versus  time  for  a 2-layer  (90/0)  laminate 
under  suddenly  applied  sinusoidal  loading 
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Figure  4.9.  Transverse  normal  stress  versus  time  for  a 2-layer  (90/0)  laminate 
under  suddenly  applied  sinusoidal  loading 
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In  Figures  4.6  and  4.7,  it  is  observed  that  the  normalized 

normal  stress  a„  at  center  top  surface  of  the  laminate  is 

close  to  the  normalized  normal  stress  5y  at  center  bottom 

surface  of  the  laminate,  except  ov  is  in  tension  and  a„  is 

x y 

in  compression.  As  shown  in  Figure  4.8,  the  variation  of 
normalized  shear  stress  xxz  is  similar  to  XyZ.  From  the 
plots,  it  is  seen  that  the  periods  of  the  transient  response 
for  w,  ax,  Oy,  t xz , and  XyZ  are  very  closely  related.  This 
fact  agrees  with  the  results  of  Putcha  and  Reddy  [45].  The 
period  for  the  normalized  transverse  normal  stress  az(a/2, 
a/2,  h/2 ) is  much  shorter  when  compared  with  others.  A 

shorter  time  step  size  (At  = 1 microsecond)  is  employed  to 
observe  the  periodic  response  of  the  normalized  transverse 
normal  stress  (Figure  4.9). 

The  normalized  central  transverse  deflection  with 
respect  to  time  for  a three-layer  simply  supported  cross-ply 
(0/90/0)  square  laminate  under  sinusoidal  loading  is  shown 
in  Figure  4.10.  Figure  4.11  contains  the  normalized  maximum 
normal  stresses,  ax  and  a , as  a function  of  time.  The 
normalized  shear  stresses,  xxz  and  XyZ,  are  shown  in  Figure 
4.12.  The  periods  for  the  normalized  deflection  and  stresses 
are  similar.  In  Figure  4.12,  the  amplitude  of  the  response 
is  larger  for  xxz  than  for  iyZ;  it  is  because  the  bending 
stiffness  is  higher  in  the  x-direction  than  in  the  y- 
direction  for  a three-layer  (0/90/0)  laminate  with  layers  of 
equal  thickness.  The  normalized  central  normal  stress 
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Figure  4.10.  Center  deflection  versus  time  for  a 3-layer  cross-ply  (0/90/0) 
laminate  under  suddenly  applied  sinusoidal  loading 


(a/2, a/2, h/2) 
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Figure  4.11.  Normal  stresses  versus  time  for  a 3-layer  (0/90/0)  laminate  under 
suddenly  applied  sinusoidal  loading 


(O.a/2,0) 
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Figure  4.12.  Transverse  shear  stresses  versus  time  for  a 3-layer  (0/90/0)  laminate 
under  suddenly  applied  sinusoidal  loading 


89 


distribution,  av,  through  the  thickness  of  the  laminate  for 
time  from  20  to  80  microseconds  is  shown  in  Figure  4.13. 

In  the  present  study,  a four-node  isoparametric  plate 
element  with  48  assumed  stress  parameters  for  each  lamina  is 
used.  Fast  convergence  is  observed;  only  a 5 x 5 mesh  is 
modeled  in  a quadrant  of  the  laminate. 
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Figure  4.13.  Through- thickness  center  normal  stress  versus  z for  a 3-layer 
(0/90/0)  laminate  under  suddenly  applied  sinusoidal  loading 
for  time  from  20  to  80  microseconds 


CHAPTER  5 


ANALYSIS  OF  LAMINATED  PLATES  UNDER 
IMPACT  DYNAMIC  LOADING 

5.1  Introduction 

From  the  work  of  previous  sections  , including  both 
static  analysis  and  dynamic  analysis,  the  accuracy  and 
versatility  of  the  present  method  have  been  validated.  It 
proves  that  the  hybrid  stress  finite  element  model  is  a 
useful  tool  in  attacking  the  complicated  structure  analysis 
problems  which  cannot  be  solved  in  the  closed  form  solution. 
In  this  chapter,  a multilayered  composite  plate  subject  to 
central  small  area  impact  loading  such  as  by  a small 
projectile  will  be  discussed. 

When  a laminated  plate  is  impacted  by  a projectile, 
flexural  waves  are  produced.  Flexural  waves  cause 
displacement  and  stress  distributions  in  the  laminated  plate 
that  vary  with  time.  The  waves  are  important  because 
experiments  conducted  by  Takeda  et  al.  [75]  indicated  that 
the  flexural  wave  propagation  is  the  principal  response  mode 
of  moderately  thin  and  thick  plates  subject  to  central 
impact,  and  it  leads  to  delamination  crack  propagation. 

The  analysis  herein  is  concerned  only  with  the  time 
that  the  laminated  plate  remains  elastic,  and  that 
delamination  has  not  yet  occurred.  The  damping  in  the 
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laminated  plate  is  assumed  to  be  negligible.  The  stresses 
and  deflections  are  presented  to  show  the  transient  response 
of  the  laminated  plate  when  impacted  by  a projectile. 
Knowledge  of  where  and  when  the  stresses  are  maximum  should 
be  helpful  in  determining  where,  when,  and  how  the 
delamination  failure  occurs.  With  the  information  of  the 
laminated  plate  behavior , the  structure  members  can  be 
constructed  to  meet  the  specific  design  requirements. 

5.2  Impact  Mechanism  and  Laminate  Characteristics 

The  impact  experimental  work  conducted  at  University  of 
Florida  has  been  described  in  chapter  1.  In  this  study,  a 
blunt-ended  steel  circular  cylinder  with  2.54  cm  in  length 
and  0.9525  cm  in  diameter  is  used.  The  projectile  has  mass 
0.014175  kg.  The  square  laminate  is  clamped  around  the  four 
edges  with  a side  length  of  14  cm.  These  dimensions  and 
projectile  mass  are  the  same  as  those  used  in  the 
experimental  works  of  Takeda  [56]. 

Two  kinds  of  material  are  used  in  the  present 
investigation.  One  is  the  Scotchply  1002  E-glass  epoxy  with 
the  following  material  properties  for  each  lamina, 

El  = 40.0  GPa  (5800000  psi) 

Et  = 8.27  GPa  (1200000  psi) 

GLT  = 4,13  GPa  (600000  Psi) 

Gtt  = 3.03  GPa  (440000  psi)  5.1 


0.26 
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p = 1.90  g/cm^ 

Tg  = 20  MPa  (2900  psi) 

where  p is  the  density  of  the  lamina,  with  units  of  mass 
per  unit  volume.  T£  is  the  transverse  tensile  strength  and 
is  assumed  to  be  equal  to  the  matrix  (epoxy)  tensile 
strength.  The  other  is  graphite  epoxy  AS-3501-6  with  the 
material  properties  as  follows: 


el  = 

142 

.73 

GPa 

(20700000  ps 

E,p  = 

13. 

79 

GPa 

( 2000000 

psi  ) 

glt 

= 4. 

64 

GPa 

( 673000 

psi ) 

Grprp 

= 4. 

14 

GPa 

(600000 

psi  ) 

VLT 

= 0. 

3 

V ipip 

= 0. 

28 

p = 1.61  g/cm 

Te  = 54  MPa  ( 7840  psi) 

In  the  present  study,  since  the  laminate  geometry, 
loading,  and  boundary  conditions  are  symmetric  with  respect 
to  axes  through  the  center  of  the  laminate,  it  is  sufficient 
to  analyze  only  one  quarter  of  the  laminate.  The  clamped 
boundary  conditions  for  the  present  model  are  shown  in 
Figure  5.1.  The  discretization  of  the  domain  of  interest  is 
arbitrary.  A typical  finite  element  mesh  is  shown  in  Figure 
5.2,  where  the  shaded  area  is  the  impact  area. 
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Figure  5.1.  Clamped  boundary  conditions  for  a 
quadrant  of  the  laminate 


Figure  5.2.  A typical  mesh  for  the  finite  element  model 
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5.3  Force  Distribution  on  a Laminated  Plate 
In  finite  element  analysis,  equivalent  nodal  forces  are 
required  in  the  computation.  However,  when  a laminate  is 
impacted  by  a projectile  only  the  velocity  of  the  projectile 
can  be  observed.  The  actual  force  distribution  between  the 
projectile  and  the  laminate  during  impact  is  unknown.  In 
this  section  methods  to  convert  the  projectile  velocity  to 
equivalent  nodal  force  distributions  will  be  discussed.  Two 
types  of  force-velocity  relation  are  simulated  in  this 
study . 

One  is  based  on  the  principle  of  linear  impulse  and 
momentum  with  assumed  contact  time  in  the  form  of 

m v = F t 5.3 

where  m is  the  mass  of  the  projectile,  v is  the  velocity  of 
the  projectile  before  impact,  and  t is  the  contact  time.  For 
a given  mass  and  velocity  of  the  projectile,  and  contact 
time,  the  time  average  of  the  time  dependent  force  F can  be 
determined  from  Equation  5.3.  Two  forms  of  the  force-time 
relation  are  modeled,  which  are 

(1)  Suddenly  applied  uniform  pulse  loading  along  the 
contact  interval  as  shown  in  Figure  5.3a. 

(2)  Suddenly  applied  triangularly  distributed  pulse 
loading  over  the  contact  interval  as  shown  in 
Figure  5.3b.  The  force  increases  linearly  rapidly 
in  a short  time  interval , then  decreases  linearly 
more  gradually  to  zero. 
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Applied  force 


F 


t 


Contact  time 


(a)  Uniformly  distributed 


(b)  Triangularly  distributed 


Figure  5.3.  Force-time  relation  based  on  the  principle 
of  linear  impulse  and  momentum 
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The  applied  force  in  case  (1)  will  not  be  equal  to  applied 
force  in  case  (2)  at  a specific  time,  but  the  overall 
impulses  for  both  cases  are  the  same,  which  leads  to  the 
highest  force  at  case  (2)  being  twice  the  uniformly 
distributed  force  in  case  (1).  At  each  time  step,  the 
loading  in  the  finite  element  analysis  is  assumed  to  be 
uniformly  distributed  over  the  circular  impact  area. 

The  other  type  of  impact  is  an  elastic  impact  including 
indentation  of  the  laminate  surface  by  the  projectile.  The 
force  between  the  projectile  and  the  laminate  during  impact 
is  assumed  to  be  governed  by  the  Hertzian  theory  in  the  form 

F = H(r  - w)P  5.4 

where  F is  the  applied  force,  H is  a material  constant,  r is 
the  displacement  of  the  projectile,  w is  the  displacement  at 
the  point  of  impact,  i.e.  the  displacement  of  the  center  of 
the  laminate,  and  p is  an  exponent.  Both  displacements  are 
measured  from  the  surface  of  the  undeformed  laminate.  The 
Hertzian  theory  of  impact  was  discussed  by  Love  [76]  and 
Goldsmith  [50].  Sun,  Sankar,  and  Tan  [58]  proposed  a simple 
static  indentation  test  to  establish  the  contact  law  for  the 
round  steel  ball  and  laminates.  They  found  that  a Hertzian 
contact  law  with  p = 1.5  is  a good  approximation  for  the 
loading-indentation  relationship.  The  Hertzian  contact  law 
was  used  successfully  by  Petersen  [77]  in  the  investigation 
of  transient  response  of  laminated  plates  subject  to  dynamic 
impact  loading.  In  his  study,  p was  set  equal  to  1.5  and 
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the  corresponding  value  of  H was  set  equal  to  100000000 
N/m  , which  were  found  to  be  appropriate.  The  Hertzian 
contact  law  with  the  same  values  of  p and  H as  Petersen's 
[77]  will  be  used  in  the  present  work.  The  contact  behavior 
is  mathematically  modeled  and  incorporated  into  the  hybrid 
stress  finite  element  program  to  perform  the  dynamic 
analysis  of  a square  laminated  plate  subject  to  impact  of  a 
steel  circular  cylinder. 

5.4  E-Glass  Epoxy  Laminates  Subject  to  Impact  Loading 

In  this  section,  the  deflections  and  stresses  in  a 
three-layer  cross-ply  (0/90/0)  laminate  subject  to  impact 
loading  will  be  examined.  The  laminate  is  made  of  E-glass 
epoxy  laminae  with  material  properties  as  given  in  Equation 
5.1.  Equal  thickness  is  assumed  in  each  layer,  the  thickness 
of  the  laminate  is  3.81  mm.  The  laminate  is  impacted  by  a 
projectile  with  mass  0.014175  kg  and  velocity  22.6  m/second. 
The  laminate  and  projectile  dimensions  are  as  described  in 
section  5.2.  The  boundary  conditions  of  the  finite  element 
model  are  given  in  Figure  5.1.  The  applied  loading  is  based 
on  Equation  5.3  with  assumed  contact  time  equal  to  300 
microseconds . 

The  transverse  deflections  of  the  center  of  the 
laminate  as  a function  of  time  are  shown  in  Figure  5.4.  The 
dotted  line  is  the  uniform  pulse  loading  along  the  contact 
interval,  and  the  solid  line  is  the  triangularly  distributed 
pulse  loading  over  the  contact  interval.  The  force-time 
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relations  are  also  included  in  Figure  5.4.  The  period  of  the 
transient  response  is  slightly  larger  for  uniformly 
distributed  loading  than  for  the  triangularly  distributed 
loading.  Regardless  the  loading  distribution  over  the 
contact  interval,  it  is  observed  that  the  central  transverse 
deflections  are  very  closely  related  each  other  as  long  as 
the  impulses  are  equal. 

The  normal  stress,  o . on  the  top  surface  of  the 
laminate,  and  shear  stress,  r xz > on  the  mid-surface  of  the 
laminate,  along  the  x axis  in  the  three-layer  laminate  under 
uniformly  distributed  loading  are  shown  in  Figures  5.5  and 
5.6  for  time  from  20  to  80  microseconds  after  impact.  The 
stress  distributions  in  the  same  laminate  under  the 
triangularly  distributed  loading  are  given  in  Figures  5.7 
and  5.8.  The  abscissa  gives  the  distance  from  the  center  of 
the  laminate.  The  center  of  the  laminate  is  set  equal  to  0 
mm,  and  the  edge  of  the  laminate  is  set  equal  to  70  mm.  The 
impact  occurs  between  0 mm  and  4.7625  mm.  Since  the  stresses 
are  assumed  within  each  element,  the  stress  continuity  at 
the  boundary  of  adjacent  elements  cannot  be  guaranteed, 
therefore  the  results  at  Gauss  points  are  chosen  to  simplify 
the  plots. 

From  Figures  5.5  and  5.7,  the  maximum  value  of  the 
normal  stress  occurs  at  the  top  center  of  the  laminate  and 
is  a compressive  stress.  The  position  of  the  maximum  tensile 
normal  stress  moves  from  the  center  towards  the  edge  of  the 
laminate  with  increasing  time.  The  motion  of  the  point  of 
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Figure  5.4.  Transverse  deflections  of  the  center  versus  time  for  a 3-layer  (0/90/0) 
E-glass  laminate  under  suddenly  applied  impact  loading 
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Figure  5.5.  Normal  stress  along  the  x axis  in  a 3-layer  E-glass  laminate  under 
uniformly  distributed  loading  for  time  from  20  to  80  microseconds 
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Figure  5.6.  Shear  stress  along  the  x axis  in  a 3-layer  E-glass  laminate 
uniformly  distributed  loading  for  time  from  20  to  80  micros 
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Figure  5.7.  Normal  stress  along  the  x axis  in  a 3-layer  E-glass  laminate 
triangularly  distributed  loading  for  time  from  20  to  80  micr 
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Figure  5.8.  Shear  stress  along  the  x axis  in  a 3-layer  E-glass  laminate 
triangularly  distributed  loading  for  time  from  20  to  80  mic 
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maximum  tensile  stress  is  caused  by  a flexural  wave 
generated  by  the  initial  impact.  The  flexural  wave  begins  at 
the  center  of  the  laminate  and  moves  outwardly.  The  shear 
stress  along  the  x axis,  Figure  5.6  or  5.8,  was  evaluated  at 
the  mid-surface  of  the  laminate.  The  maximum  value  of  the 
shear  stress  occurs  near  the  outer  edge  of  the  impact  area. 
The  shear  stress  at  the  center  of  the  laminate  is  equal  to 
zero  because  of  biaxial  symmetry.  The  wave  motion  of  the 
shear  stress  distribution  which  is  caused  by  a flexural  wave 
propagating  from  the  area  of  impact  toward  the  edge  of  the 
laminate  is  also  observed. 

As  previously  mentioned,  the  central  transverse 
deflection  of  a 3-layer  cross-ply  laminate  under  uniformly 
distributed  loading  was  in  good  agreement  with  those  of 
triangularly  distributed  loading.  However,  the  stresses  are 
significantly  dependent  on  the  history  of  the  applied 
loading.  For  a uniformly  distributed  loading,  Figure  5.6, 
the  maximum  value  of  the  shear  stress  is  almost  the  same 
with  increasing  time.  In  the  case  of  triangularly 
distributed  loading,  Figure  5.8,  the  maximum  value  of  the 
shear  stress  is  increasing  gradually  corresponding  to  the 
increase  of  the  applied  loading.  Similar  behavior  is  also 
noticed  in  the  normal  stress  distributions. 

In  the  present  study,  the  central  transverse  deflection 
w increases  as  time  t increases  until  t reaches  300 
microseconds  as  shown  in  Figure  5.4.  The  stresses  however, 
during  this  time  interval  do  not  in  general  increase  as  the 
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transverse  deflection.  These  facts  can  be  observed  in 
Figures  5.9  and  5.10  which  show  the  normal  stress  and 
transverse  shear  stress  along  the  x axis  of  a 3-layer 
laminate  under  triangularly  distributed  loading  for  time 
from  100  to  400  microseconds. 

The  deflection  and  stresses  in  the  previous  example 
were  calculated  with  the  assumptions  of  contact  time  and 
loading  distribution.  In  the  present  example  an  elastic 
impact  is  simulated.  The  elastic  impact  allows  the 
projectile  to  indent  the  laminate  and  rebound.  The  applied 
force  between  the  projectile  and  the  laminate  during  impact 
is  assumed  to  be  governed  by  the  Hertzian  impact  law 
(Equation  5.4).  p is  set  equal  to  1.5  and  H is  chosen  to  be 
100000000  N/m^'^,  as  discussed  in  section  5.3. 

Figure  5.11  contains  plots  of  the  central  transverse 
deflection  of  a 3-layer  (0/90/0)  E-glass  epoxy  laminate  as  a 
function  of  time  determined  experimentally  and  by  the  hybrid 
stress  finite  element  method.  The  laminate  thickness  is  4.29 
mm  and  the  projectile  velocity  is  equal  to  22.6  m/second. 
The  period  of  response  from  calculation  is  shorter  for  the 
case  of  140  mm  x 140  mm  square  laminate  than  the 
experimentally  determined  result.  It  is  also  noticed  that 
the  deflection  evaluated  by  the  hybrid  stress  finite  element 
method  with  the  side  length  equal  to  140  mm  agrees  with  the 
deflection  determined  experimentally  for  time  interval  less 
than  200  microseconds.  As  time  increases,  the  computed  and 
experimental  deflections  do  not  agree  closely  because  that 
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Figure  5.9.  Normal  stress  along  the  x axis  in  a 3-layer  E-glass  laminate  under 

triangularly  distributed  loading  for  time  from  100  to  400  microseconds 
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Figure  5.10.  Shear  stress  along  the  x axis  in  a 3-layer  E-glass  laminate  under 

triangularly  distributed  loading  for  time  from  100  to  400  microseconds 
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the  delamination  is  initiated  at  the  interface  of  the 
laminate  and  the  experimental  boundary  condition  is  not 
truly  clamped  around  the  four  edges. 

The  metal  frame  which  holds  the  impacted  laminate  will 
allow  the  laminate  to  move  in  the  horizontal  direction  and 
rotate  a small  amount.  These  effects  cause  that  the 
effective  side  length  of  the  laminated  plate  is  larger  than 
the  distance  between  the  supports.  With  the  same  thickness 
of  the  laminate  and  the  projectile  velocity,  the  transverse 
deflection  evaluated  with  side  length  equal  to  170  mm  is 
also  presented.  A better  agreement  is  achieved  between  the 
computed  and  experimental  deflections. 

In  the  following  investigations,  the  effort  will  be 
concentrated  on  the  numerical  analysis.  A perfect  clamped 
boundary  condition  is  assumed  around  the  edges  and  the  side 
length  will  be  set  equal  to  140  mm  in  the  computation.  The 
central  transverse  deflection  of  a 3-layer  E-glass  epoxy 
laminate  as  a function  of  time  with  Hertzian  impact  is  shown 
in  Figure  5.12.  The  laminate  has  the  same  material 
properties  and  geometry  ( thickness=3 . 8 1 mm)  as  the  previous 
laminate  used  in  the  analysis  with  the  assumed  contact  time 
and  applied  force  distribution.  The  projectile  velocity  is 
equal  to  22.6  m/second.  The  results  of  previous 
investigation.  Figure  5.4,  are  also  included  for  comparison. 
The  differences  among  them  for  time  less  than  600 
microseconds  were  caused  by  the  assumptions  of  the  magnitude 
of  the  applied  force  and  the  contact  time.  For  time  greater 


Applied  force-time  distribution 
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Figure  5.12.  Transverse  deflections  of  the  center  versus  time  for  a 3-layer  E-glass 
laminate  under  Hertzian  impact,  h=3.81  mm,  v=22.6  m/second 
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than  600  microseconds,  the  deflection  computed  with  Hertzian 
impact  is  smaller  than  those  of  the  previous  example.  This 
is  owing  to  the  fact  that  the  Hertzian  impact  law  allows 
that  the  projectile  recontacts  the  laminate  and  force  occurs 
again  between  them.  This  effect  can  be  observed  from  the 
calculated  force  between  projectile  and  the  laminate  with 
Hertzian  impact  law,  which  is  also  given  in  Figure  5.12.  The 
first  impact  occurs  from  0 to  approximately  260 
microseconds,  at  which  time,  the  projectile  loses  contact 
with  the  laminate.  The  projectile  recontacts  the  laminate  at 
approximately  530  microseconds.  These  observations  agree 
with  the  results  of  Takeda 's  experimental  works  [56].  For  a 
higher  velocity  of  the  projectile,  a longer  contact  time  is 
observed.  Because  the  difficulties  exist  in  determining  the 
contact  time  and  describing  the  shape  of  the  applied 
loading,  it  is  believed  that  incorporating  the  Hertzian 
impact  law  into  the  finite  element  program  to  analyze  the 
impact  problems  is  appropriate  and  will  give  more  accurate 
laminate  response.  The  Hertzian  impact  law  will  be  applied 
in  the  following  examples. 

The  transverse  deflections  along  the  x axis  and  the 
y axis  of  the  same  laminate  under  Hertzian  impact  for  time 
from  40  to  200  microseconds  are  shown  in  Figures  5.13  and 
5.14.  The  results  show  the  flexural  wave  motion  along  the  x 
axis  and  the  y axis  from  the  center  to  the  edge  of  the 
laminate.  Since  EL  is  greater  than  ET,  the  bending  stiffness 
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22  for  a 3-layer  cross-ply  (0/90/0) 


is  greater  than  D 
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Figure  5.13.  Transverse  deflection  along  the  x axis  in  a 3-layer  E-glass  laminate 
under  Hertzian  impact  for  time  from  40  to  200  microseconds, 
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laminate.  This  difference  in  bending  stiffness  causes  the 
flexural  wave  to  move  faster  in  the  x direction  than  in  the 
y direction.  For  example,  observe  the  deflection  at  time 
equal  to  80  microseconds,  the  distance  from  the  center  to 
the  nearest  point  of  zero  deflection  is  36  mm  in  the  x 
direction  and  28  mm  in  the  y direction.  At  a particular 
time,  the  transverse  deflection  is  larger  at  a specified 
point  along  the  x axis  than  at  a point  at  the  same  distance 
along  the  y axis.  The  difference  between  them  is  owing  to 
the  anisotropy  exhibited  by  a 3-layer  cross-ply  laminate. 

The  normal  stress,  a , transverse  normal  stress,  az , 
and  transverse  shear  stress,  x , along  the  x axis  in  the  3- 
layer  laminate  are  shown  in  Figures  5.15,  5.16,  and  5.17  for 
time  from  20  to  80  microseconds.  The  normal  stress,  a , 
which  is  evaluated  at  the  upper  surface  of  the  90°  layer 
(2/3  laminate  thickness  above  the  bottom  of  the  laminate) 
and  the  transverse  shear  stress,  TyZr  along  the  y axis  for 
time  from  20  to  80  microseconds  are  shown  in  Figures  5.18 
and  5.19.  As  the  observation  for  the  transverse  deflection, 
in  which  the  flexural  wave  moves  faster  in  the  x direction 
than  in  the  y direction,  the  stress  wave  is  also  noticed 
moving  in  a similar  fashion.  The  maximum  value  of  the  normal 
stress,  ax,  along  the  x axis  is  greater  than  the  maximum 
value  of  the  normal  stress,  a , along  the  y axis  because  the 
modulus  of  elasticity  is  greater  in  the  x direction  than  in 
the  y direction.  The  maximum  value  of  the  shear  stress,  Txz / 
at  the  outer  edge  of  the  impact  zone  on  the  x axis  is 
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Figure  5.15.  Normal  stress  along  the  x axis  in  a 3-layer  E-glass  laminate  under 

Hertzian  impact  for  time  from  20  to  80  microseconds,  v=22.6  m/second 
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Figure  5.16.  Transverse  normal  stress  along  the  x axis  in  a 3-layer  E-glass  laminate 

under  Hertzian  impact  for  time  from  20  to  80  microseconds,  v=22.6  m/second 
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Figure  5.18.  Normal  stress  along  the  y axis  in  a 3-layer  E-glass  laminate  under 

Hertzian  impact  for  time  from  20  to  80  microseconds,  v=22.6  m/second 
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Figure  5.19.  Shear  stress  along  the  y axis  in  a 3-layer  E-glass  laminate  under 

Hertzian  impact  for  time  from  20  to  80  microseconds,  v=22.6  m/second 
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slightly  higher  than  the  maximum  value  of  the  shear  stress, 
TyZ,  on  the  y axis.  The  transverse  normal  stress  (Figure 
5.16),  az , which  is  evaluated  at  the  top  surface  of  the 
laminate  is  significantly  dependent  on  the  applied  force 
computed  from  the  Hertzian  impact  law  (Figure  5.12).  Outside 
the  impact  region,  the  transverse  normal  stress  approaches 
to  zero.  Therefore,  transverse  normal  stress  on  the  outside 
of  the  impact  region  will  not  have  a significant  effect  on 
the  initiation  of  delamination  of  laminate  when  it  is 
subject  to  impact  loading. 

Figures  5.20  and  5.21  contain  the  normal  stress,  a„, 
and  transverse  shear  stress,  Txz>  along  the  x axis  for  time 
from  100  to  400  microseconds.  It  is  observed  that  the 
transverse  shear  stress  near  the  impact  region  is  much 
greater  for  time  equal  to  100  microseconds  than  for  time 
equal  to  300  microseconds.  This  fact  is  because  that  at  300 
microseconds  the  projectile  has  lost  contact  with  the 
laminate  and  no  force  occurs  between  them  as  shown  in 
Figure  5.12.  It  is  clear  that  the  transverse  shear  stress  is 
dependent  on  the  history  of  the  applied  loading. 

Figure  5.22  contains  plots  of  the  central  transverse 
deflection  versus  time  for  a 3-layer  cross-ply  (0/90/0) 
laminate  under  three  different  projectile  velocities,  20 
m/second,  30  m/second,  and  40  m/second,  respectively.  The 
amplitude  of  the  central  transverse  deflection  increases 
with  increasing  projectile  velocity  and  is  almost 
proportional  to  the  projectile  velocity,  but  the  period  of 
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Figure  5.20.  Normal  stress  along  the  x axis  in  a 3-layer  E-glass  laminate  under 

Hertzian  impact  for  time  from  100  to  400  microseconds,  v=22.6  m/second 
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Figure  5.22.  Transverse  deflections  of  the  center  versus  time  for  a 3-layer  E-glass 
laminate  under  Hertzian  impact  with  three  projectile  velocities, 
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the  dynamic  response  remains  the  same.  It  is  seen  that  the 
projectile  velocity  does  not  have  an  effect  on  the  period  of 
laminate  response.  The  normal  stress,  crx,  and  transverse 
shear  stress,  Txz>  along  the  x axis  in  the  3-layer  laminate 
with  projectile  velocity  equal  to  40  m/second  are  shown  in 
Figures  5.23  and  5.24  for  time  from  20  to  80  microseconds. 
As  expected,  stress  increases  with  increasing  projectile 
velocity.  At  a particular  time,  the  distance  from  the  center 
to  the  nearest  point  of  the  laminate  where  the  stress  is 
zero  is  the  same  as  in  the  case  for  projectile  velocity 
equal  to  22.6  m/second  (Figures  5.15  and  5.17).  It  is  clear 
that  flexural  wave  velocities  are  little  affected  by  the 
projectile  velocity  for  time  less  than  100  microseconds 
after  impact. 

5.5  Graphite  Epoxy  Laminates  Subject  to  Impact  Loading 
In  the  previous  section,  the  deflections  and  stresses 
in  a 3-layer  cross-ply  E-glass  epoxy  laminate  subject  to 
different  loadings  and  projectile  velocities  were  presented. 
In  this  section,  the  dynamic  response  of  graphite  epoxy 
laminates  will  be  discussed.  The  laminate  has  the  same 
dimensions  as  the  E-glass  epoxy  laminate.  The  material 
properties  of  graphite  epoxy  were  given  in  Equation  5.2.  In 
all  the  examples,  the  applied  force  between  the  projectile 
and  the  laminate  is  modeled  by  the  Hertzian  impact  law. 
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Figure  5.23.  Normal  stress  along  the  x axis  in  a 3-layer  E-glass  laminate  under 
Hertzian  impact  for  time  from  20  to  80  microseconds,  v=40  m/second 
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Figure  5.25  shows  the  maximum  central  transverse 
deflection  as  a function  of  time  for  a symmetric  three-layer 
(0/90/0)  square  laminate  with  projectile  velocity  equal  to 
22.6  m/second.  Because  of  the  greater  stiffness  the 
amplitude  of  the  central  deflection  is  smaller  and  the 
period  of  response  is  shorter  in  the  graphite  epoxy  than  in 
the  E-glass  epoxy  (Figure  5.12)  for  the  same  laminate 
geometry  and  projectile  velocity.  It  is  also  noted  that 
higher  bending  stiffness  in  the  x direction  than  in  the  y 
direction  causes  the  flexural  waves  to  move  faster  in  the  x 
direction  than  in  the  y direction.  This  effect  can  be  seen 
in  Figures  5.26  and  5.27  which  describe  the  transverse 
deflection  along  the  x axis  and  the  y axis  for  time  from  40 
to  200  microseconds  after  impact. 

The  normal  stress,  crx,  and  transverse  shear  stress, 
t , along  the  x axis  in  the  3-layer  graphite  laminate  with 
projectile  velocity  equal  to  22.6  m/second  are  shown  in 
Figures  5.28  and  5.29  for  time  from  20  to  80  microseconds 
after  impact  and  in  Figures  5.30  and  5.31  for  time  from  100 
to  400  microseconds.  The  maximum  normal  stress,  Oy,  and 
transverse  shear  stress,  tyZ,  along  the  y axis  in  the  3- 
layer  laminate  for  the  same  impact  are  shown  in  Figures  5.32 
and  5.33  for  time  from  20  to  80  microseconds.  As  in  the  case 
of  3-layer  E-glass  epoxy  laminate,  the  same  characteristics 
in  the  stress  distributions  are  observed  in  the  present 
example,  except  the  stress  magnitude  and  the  flexural  wave 
velocity . 
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Figure  5.25.  Transverse  deflection  of  the  center  versus  time  for  a 3-layer  graphite 
laminate  under  Hertzian  impact,  h=3.81  mm,  v=22.6  m/second 
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Figure  5.27.  Transverse  deflection  along  the  y axis  in  a 3-layer  graphite  laminate 
under  Hertzian  impact  for  time  from  40  to  200  microseconds, 
v=22.6  m/second 


Distance  from  impact  (mm) 


132 


Figure  5.28.  Normal  stress  along  the  x axis  in  a 3-layer  graphite  laminate 
Hertzian  impact  for  time  from  20  to  80  microseconds,  v=22.6  m 
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Figure  5.29.  Shear  stress  along  the  x axis  in  a 3-layer  graphite  laminate  under 

Hertzian  impact  for  time  from  20  to  80  microseconds,  v=22.6  m/second 
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Figure  5.30.  Normal  stress  along  the  x axis  in  a 3-layer  graphite  laminate  under 

Hertzian  impact  for  time  from  100  to  400  microseconds,  v=22.6  m/second 
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Figure  5.32.  Normal  stress  along  the  y axis  in  a 3-layer  graphite  laminate 
Hertzian  impact  for  time  from  20  to  80  microseconds,  v=22.6  m 


137 


o o c o 
cn  «a-  m co 


o 

•* 

>> 

* 

o 

*3 

N 

Q- 

>> 

s: 

H 


<X3 

G 


g 

0 

o 

O 

•a 

0) 

c 

in 

G \ 

6 

aj 

■p 

kO 

(0 

• 

G cn 

•H 

IN 

e 

II 

ra 

> 

l"H 

*» 

a) 

in 

-P  T3 

•H 

G 

g 

0 

& U 

cd 

(1) 

G 

in 

CP  0 

G 

G 

o 

<U 

•H 

>1  a 

id 

i — i 

o 

1 

CO 

m 

0 

id 

-p 

C 

o 

•rH 

CN 

in 

a 

•H 

0 

X 

G 

(0 

>P 

>1  0) 

a 

a> 

■rH 

G 

-P 

P 

G 

Cn  0 

G 

fal 

0 

r— 1 

■p 

Id 

O 

id 

in 

a 

in 

a 

<u 

•rH 

G 

+j 

c 

in 

id 

*H 

G 

N 

id 

-P 

QJ 

G 

G 

0) 

to 

« 

n 

n 

• 

in 

<u 

g 

3 

Cn 

■H 

fa 


138 


The  effect  of  laminate  thickness  on  the  central 
transverse  deflection  is  shown  in  Figure  5.34  for  a 3-layer 
cross-ply  (0/90/0)  square  graphite  laminate  under  Hertzian 
impact.  Three  various  laminate  thickness  are  modeled  in  the 
present  investigation,  3.81  mm,  5.715  mm,  and  7.62  mm, 
respectively.  The  projectile  velocity  is  equal  to  40 
m/second.  It  is  observed  that  the  amplitude  and  the  period 
of  response  of  the  central  transverse  deflection  are 
decreasing  with  increasing  laminate  thickness.  The  thickness 
of  the  laminate  has  a significant  effect  on  the  transverse 
deflection  and  normal  stress  at  the  top  center  of  the 
laminate.  For  a thicker  laminate,  the  transverse  deflection 
and  normal  stress,  ax,  reduce  a great  amount  but  the 
transverse  shear  stress,  txz , is  only  slightly  less  than  the 
thinner  laminate.  It  is  also  noticed  that  the  flexural  wave 
velocity  is  affected  by  the  thickness  of  the  laminate.  The 
wave  motion  of  stress  distribution  is  faster  in  the  thicker 
laminate  than  in  the  thinner  laminate. 

The  final  example  considered  is  the  effect  of  laminate 
anisotropy  on  the  dynamic  response  of  cross-ply  laminate. 
Three  types  of  laminations,  3-layer  cross-ply  (0/90/0),  4- 
layer  cross-ply  (0/90/0/90),  and  5-layer  cross-ply 
(0/90/0/90/0),  are  used  in  the  analysis.  Each  laminate  has 
the  same  thickness  of  3.81  mm  and  constructed  with  layers  of 
equal  thickness.  The  central  transverse  deflection  versus 
time  for  the  three  laminates  with  the  same  projectile 
velocity  of  22.6  m/second  are  shown  in  Figure  5.35.  Since 


h=3 . 81  mm 
h=5 . 71 5 mi 
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Figure  5.34.  Transverse  deflections  of  the  center  versus  time  for  a 3-layer  graphite 
laminate  under  Hertzian  impact  with  three  laminate  thickness, 
v=40  m/second 
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greater  anisotropy  is  exhibited  in  3-layer  and  5-layer 
symmetric  laminates,  the  amplitude  of  the  transverse 
deflection  is  smaller  when  compared  with  the  4-layer  anti- 
symmetric laminate.  The  normal  stress,  o and  transverse 
shear  stress,  Txz,  along  the  x axis  for  the  5-layer  laminate 
are  shown  in  Figures  5.36  and  5.37  for  time  from  40  to  200 
microseconds.  The  bending  stiffness  D-^  is  lower  in  the  5- 
layer  laminate  than  in  the  3-layer  laminate.  The  flexural 
wave  motion  which  can  be  observed  in  the  stress 
distributions  is  slightly  slower  in  the  5-layer  laminate 
than  in  the  3-layer  laminate  (Figures  5.28  and  5.29). 

5.6  Summary  Remarks 

The  deflections  and  stresses  of  laminated  fiber- 
reinforced  composites,  E-glass  epoxy  and  graphite  epoxy, 
have  been  analyzed  and  presented  graphically.  Two  types  of 
force  distribution  were  modeled  in  the  hybrid  stress  finite 
element  analysis;  one  is  based  on  the  principle  of  linear 
impulse  and  momentum  with  assumed  contact  time  and  the  other 
is  based  on  the  Hertzian  impact  law.  For  the  same  contact 
time  and  momentum  of  impact,  the  force  distribution  over  the 
contact  interval  had  little  effect  on  the  central  transverse 
deflection,  but  had  a significant  influence  on  the  stress 
distribution.  In  general,  the  contact  time  and  the 
distribution  of  applied  force  were  not  known  for  an  impact 
problem;  therefore  an  appropriate  contact  law  called 
Hertzian  impact  was  incorporated  into  the  finite  element 
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Figure  5.35.  The  effect  of  laminate  anisotropy  on  the  central  transverse  deflection 
versus  time  for  3-layer,  4-layer,  and  5-layer  graphite  laminate  under 
Hertzian  impact,  v=22.6  m/second 


Distance  from  impact  (mm) 
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Figure  5.36.  Normal  stress  along  the  x axis  in  a 5-layer  graphite  laminate 
Hertzian  impact  for  time  from  40  to  200  microseconds,  v=22.6 
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Figure  5.37.  Shear  stress  along  the  x axis  in  a 5-layer  graphite  laminate  under 

Hertzian  impact  for  time  from  40  to  200  microseconds,  v=22.6  m/second 
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program  to  perform  the  analysis.  The  force  distribution 
determined  by  the  Hertzian  impact  showed  that  the  projectile 
impacted  the  laminated  plate,  lost  contact,  and  then 
recontacted,  which  agreed  with  the  experimental 
observations.  It  is  believed  that  the  Hertzian  impact  law 
can  describe  the  contact  behavior  better  and  give  more 
accurate  solutions  than  the  assumed  contact  time. 

The  stress  distributions  along  the  x axis  and  the  y 
axis  showed  evidence  of  flexural  wave  motion  produced  by 
impact  on  the  laminated  plate.  The  flexural  wave  started 
from  the  center  of  impact  and  propagated  outwardly.  For  an 
anisotropic  laminated  plate,  the  flexural  wave  propagated 
faster  in  the  direction  in  which  the  bending  stiffness  is 
the  greatest.  Higher  stresses  were  also  exhibited  in  the 
direction  having  higher  stiffness.  The  transverse  normal 
stress  and  transverse  shear  stresses  were  dependent  on  the 
history  of  the  applied  loading.  After  the  time  when  the 
projectile  separated  from  the  laminated  plate,  these 
stresses  decreased  significantly. 

The  deflections  and  stresses  were  proportional  to  the 
projectile  velocity,  although  the  period  of  response 
remained  the  same.  The  projectile  velocity  had  little  effect 
on  the  velocity  of  the  flexural  wave  for  time  less  than  100 
microseconds.  For  a thicker  laminated  plate,  the  amplitude 
of  the  central  transverse  deflection  and  the  period  of 
response  were  smaller.  Laminate  thickness  also  had  an 
influence  on  the  flexural  wave  velocity.  The  flexural  wave 
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moved  toward  the  edge  of  laminated  plate  faster  in  the 
thicker  laminate  than  in  the  thinner  laminate.  The  magnitude 
of  stresses  were  also  affected  by  the  thickness  of  the 
laminate.  The  transverse  deflection  was  affected  by  the 
laminate  anisotropy.  The  amplitude  of  the  central  transverse 
deflection  was  smaller  in  the  3-layer  and  5-layer  symmetric 
laminates  which  possess  greater  anisotropy  than  in  the  4- 
layer  anti-symmetric  laminate.  The  flexural  wave  propagated 
faster  in  the  laminate  which  has  higher  bending  stiffness. 

The  deflection  and  stress  analysis  for  the  present 
impact  problems  was  based  on  the  assumption  of  linear 
elasticity.  It  was  assumed  that  no  failures  will  occur 
within  the  material.  However,  experimental  studies  of  impact 
conducted  by  Takeda  [56]  showed  that  delamination  usually 
exists  at  the  interface  between  layers.  In  1982,  Takeda  et 
al.  [75]  showed  that  impact  causes  the  delamination  crack  to 
propagate.  They  observed  a delamination  crack  beginning  near 
the  point  of  impact  and  propagating  outwardly.  It  is  seen 
that  the  flexural  wave  which  causes  the  stresses  to 
propagate  may  be  the  major  reason  for  the  delamination 
propagation . 

Composite  laminates  are  typically  constructed  by 
stacking  the  individual  laminae  together  with  the  material 
called  matrix.  Since  the  matrix  has  low  strength  in 
comparison  to  the  fiber  strength,  the  interlaminar 
transverse  shear  stresses  have  more  important  influence  for 
laminated  plates  than  for  isotropic  plates  and  usually  cause 
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delamination  cracks  at  the  interfaces.  Pagano  and  Pipes  [78] 
suggested  that  interlaminar  normal  stress  is  also 
responsible  for  the  delamination  cracks. 

The  present  hybrid  stress  finite  element  model  is 
capable  of  providing  predictions  of  the  interlaminar 
stresses  which  actually  cause  delamination  cracks  at  the 
interfaces.  It  is  reasonable  to  incorporate  the  interlaminar 
stresses  into  a suitable  failure  criterion  to  perform  the 
failure  analysis  and  observe  the  delamination  crack.  The 
delamination  crack  propagation  in  a laminated  composite 
plate  is  a very  complicated  subject.  After  the  delamination 
crack  occurs,  the  stiffness  of  the  laminated  plate  has  been 
changed,  which  will  lead  to  the  topic  of  fracture  mechanics. 
The  assumptions  in  the  present  model  may  not  be  able  to 
predict  the  delamination  crack  propagation  adequately,  but 
it  can  provide  some  valuable  information  about  the 
initiation  of  the  delamination  crack  and  it  is  hoped  lead  to 
further  understanding  of  the  delamination  propagation 
mechanism. 


CHAPTER  6 


FAILURE  ANALYSIS  OF  LAMINATED  PLATES 
6.1  Introduction 

Because  of  the  many  advantages  existing  in  the 
composite  materials,  the  application  of  laminated  plates  has 
grown  rapidly  in  recent  years.  Frequently,  failures  which 
include  fiber  breaking,  matrix  cavitation,  crack  propagation 
and  delamination  crack  are  observed  in  the  composite 
structural  elements.  Therefore,  an  appropriate  failure 
criterion  is  recommended  to  analyze  the  strengths  of  the 
structural  elements  such  that  it  will  meet  the  specific 
design  requirements.  There  have  been  many  failure  criteria 
in  existence,  such  as  maximum  stress  criterion,  maximum 
strain  criterion,  Tsai-Hill  criterion  [79],  and  Tsai-Wu 
criterion  [80].  These  criteria  may  not  possibly  describe  the 
actual  fracture  mechanisms,  but  they  can  at  least  be  used  as 
tools  for  the  material  characterization. 

In  the  present  investigation,  only  the  delamination 
crack  at  the  interface  will  be  considered.  A simple 
quadratic  failure  criterion  will  be  applied  to  perform  the 
analysis.  This  criterion  only  takes  into  account  the 
interlaminar  transverse  shear  stresses  and  transverse  normal 
stress,  which  are  believed  to  be  the  major  causes  for  the 
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delamination  crack.  Since  matrix  is  the  only  material  and 
bonding  agent  between  laminae  and  there  is  no  fiber 
breaking,  the  problem  will  be  treated  as  a failure  analysis 
of  an  isotropic  material. 

The  stress  distributions  at  the  interface  for  both  E- 
glass  epoxy  and  graphite  epoxy  laminated  plates  subject  to 
impact  loading  will  be  computed  by  using  the  hybrid  stress 
finite  element  method.  The  results  are  incorporated  into  the 
failure  criterion  to  examine  the  failure  envelope  at  a 
particular  interface.  For  simplicity  and  smoothness  of  the 
results,  cubic  spline  curve  fitting  technique  [81]  will  be 
adopted  to  plot  the  results  which  will  give  when,  how,  and 
where  the  delamination  crack  occurs.  Since  the  material  was 
assumed  within  the  linear  elastic  range  in  the  stress 
calculation,  the  stiffness  of  the  laminated  plate  does  not 
decrease  after  initial  delamination  occurs.  The  results 
presented  herein  cannot  explain  the  actual  delamination 
crack  propagation  phenomenon  and  the  final  shape  of  the 
delamination  area.  They  can  only  provide  some  information 
about  the  initiation  and  the  characteristics  of  the 
delamination  mechanism. 

6.2  Quadratic  Failure  Criterion 

As  previously  mentioned,  the  delamination  crack  is  the 
main  concern  in  the  present  investigation  and  interlaminar 
transverse  shear  stresses,  xxz  and  TyZ , and  transverse 
normal  stress,  qz,  are  the  major  causes  for  the  delamination 
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crack.  The  present  assumed  criterion  states  that  in  the 
interface  of  the  laminated  plate  there  exists  a failure 
envelope  in  the  form  of 


f (T) 


1 


6.1 


where  Y and  S are  the  tensile  strength  and  pure  shear 
strength  of  the  matrix,  respectively,  where 
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+ X 
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It  was  assumed  that 


Y2  = 3S2 


6.3 


in  analogy  with  the  Von  Mises  yield  condition  for  an 
isotropic  material  [60],  although  this  may  not  be  the  best 
choice  for  an  interlaminar  brittle  failure  criterion. 
Equation  6.1  becomes 


f (T) 
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6.4 


The  failure  criterion,  Equation  6.4,  is  the 
simplification  of  Mises  criterion  [60]  if  inplane  normal 
stresses  and  shear  stress  are  ignored.  It  is  also  noticed 
that  the  criterion  can  be  derived  from  Tsai-Wu  criterion 
[80]  for  the  isotropic  material.  Equation  6.4  is  not  a 
general  theory  for  the  failure  analysis  of  laminated 
structure;  it  can  only  be  applied  to  the  delamination  crack 
analysis.  Some  numerical  examples  will  be  discussed  in  the 


next  section. 
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6.3  Delamination  Cracks  of  Laminated  Plates 

Numerical  results  obtained  by  using  the  hybrid  stress 
finite  element  method  and  the  failure  criterion  developed 
herein  are  presented  for  the  delamination  analysis  of  a 
laminated  plate  subject  to  the  impact  loading.  The  laminated 
plate  dimension  and  boundary  conditions,  and  projectile  mass 
and  dimension  are  the  same  as  those  used  in  the  previous 
chapter . 

A failure  envelope  of  a three-layer  cross-ply  (0/90/0) 
E-glass  epoxy  square  laminate  with  layers  of  equal  thickness 
subject  to  impact  loading  determined  experimentally  [56]  and 
compared  a contour  of  the  failure  criterion  calculated  by 
the  hybrid  stress  finite  element  method  are  shown  in  Figure 
6.1,  where  x and  y indicate  the  fiber  direction  and 
transverse  direction,  respectively.  The  laminate  thickness 
is  equal  to  3.81  mm  and  the  material  properties  of  E-glass 
epoxy  were  given  in  Equation  5.1  where  the  tensile  strength 
of  epoxy  was  equal  to  20  MPa.  The  projectile  velocity  is 
equal  to  75.1  m/second.  The  contact  force  is  governed  by  the 
Hertzian  impact  law.  Owing  to  the  high  projectile  velocity, 
the  experimental  result  showed  that  the  fiber  breaking, 
transverse  cracks,  and  delamination  cracks  were  exhibited  in 
the  laminated  plate  [56].  Since  in  the  present  failure 
criterion  assumptions,  the  material  stiffness  does  not 
decrease  after  the  initial  delamination  crack  occurred,  the 
computed  envelope  cannot  match  the  experimental  result. 
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Figure  6.1.  Delamination  failure  envelope  of  a 3-layer 
E-glass  laminate  determined  experimentally 
and  by  using  hybrid  stress  finite  element 
method,  h=3.81  mm,  v=75.1  m/second 
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In  the  following  examples,  low  projectile  velocities 
will  be  modeled.  Because  of  the  lack  of  experimental 
results,  only  numerical  results  are  presented.  The  same 
laminate  geometry  and  boundary  conditions,  but  with  the 
projectile  velocity  being  reduced  to  22.6  m/second,  will  be 
investigated  next.  Figure  6.2  shows  the  contours  of  function 
f value  (left-hand  side  of  Equation  6.4)  for  both  interfaces 
at  time  equals  40  microseconds.  At  that  moment  a complete 
curve  for  f=l  is  observed  at  the  interfaces.  At  times  before 
40  microseconds,  some  curves  for  f = l are  noticed,  but  they 
do  not  make  complete  curves.  The  interfaces  are  numbered 
from  bottom  to  top  of  the  laminate,  which  means  that  the 
interface  near  the  impact  side  will  be  the  second  interface. 
Because  of  the  biaxial  symmetry,  only  the  result  in  the 
first  quadrant  of  the  laminate  is  plotted.  The  result  shows 
that  the  value  of  function  f is  decreasing  with  distance 
from  the  center  of  the  laminate.  This  fact  agrees  with  the 
observations  of  transverse  shear  stress  distributions  in 
chapter  5,  which  described  that  the  transverse  shear  stress 
is  maximum  near  the  boundary  of  the  impact  area.  For  the 
value  of  function  f equal  to  unity,  it  implies  that  the 
delamination  crack  has  occurred  inside  the  region  enclosed 
by  the  curve.  The  delamination  crack  will  not  occur  outside 
the  region,  because  that  the  contours  have  the  values  of 
function  f less  than  unity.  It  is  concluded  that  the 
delamination  crack  will  initiate  near  the  impact  area  of  the 
laminate.  The  results  for  both  interfaces  are  very  similar 
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(a)  First  interface 


(b)  Second  interface 

Contours  of  failure  value  for  both  interfaces 
of  a 3-layer  E-glass  laminate  under  Hertzian 
impact  at  time  equal  to  40  microseconds, 
h=3.81  mm,  v=  2 2 .6  m/second 


Figure  6.2. 
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to  each  other.  This  indicates  that,  based  on  this  failure 
criterion,  the  delamination  cracks  for  both  interfaces  may 
possibly  initiate  at  the  same  moment. 

For  the  same  example,  Figure  6.3  contains  the  plots  of 
the  f=l  contour  propagations  for  both  interfaces  for  time 
from  40  to  100  microseconds.  The  result  shows  that  the  f = l 
contour  is  moving  outwardly  toward  the  edge  of  the  laminate 
after  impact.  Because  the  stiffness  of  the  material  does  not 
decrease  after  the  initial  delamination  occurred,  the  result 
cannot  explain  the  actual  phenomenon  of  the  delamination 
crack  propagation.  It  describes  that  the  delamination  crack 
is  possibly  to  propagate  in  a laminated  plate  when  it  is 
under  localized  impact.  Note  that  the  f = l contour 
propagations  in  both  interfaces  are  very  closely  related. 

The  f=l  contour  propagations  in  a three-layer  cross-ply 
(0/90/0)  E-glass  epoxy  laminate  under  Hertzian  impact  with 
projectile  velocityis  equal  to  30.0  m/second  is  shown  in 
Figure  6.4  for  time  from  30  to  90  microseconds.  The 
laminated  plate  geometry  and  material  properties  are  the 
same  as  the  previous  example.  As  expected,  the  greater 
interlaminar  stresses  were  exhibited  for  the  higher 
projectile  velocity,  the  delamination  crack  is  initiated 
earlier  for  the  higher  projectile  velocity  than  for  the 
lower  projectile  velocity.  The  shape  of  the  delamination 
crack  is  also  affected  by  the  projectile  velocity,  larger 
f=l  contour  propagations  are  shown  for  the  higher  projectile 
velocity . 
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(a)  First  interface 
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(b)  Second  interface 

Figure  6.3.  The  f=l  contour  propagations  for  both  interfaces 
of  a 3-laver  E-glass  laminate  under  Hertzian 
impact  for  time  from  40  to  100  microseconds, 
h=3.81  mm,  v=22.6  m/second 
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Figure  6.4.  The  f=l  contour  propagations  for  both  interfaces 
of  a 3-layer  E-glass  laminate  under  Hertzian 
impact  for  time  from  30  to  90  microseconds, 
h=3.81  mm,  v=30  m/second 
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The  delamination  analysis  for  a three-layer  cross-ply 
(0/90/0)  graphite  epoxy  laminate  will  be  presented  next.  The 
laminated  plate  has  the  same  dimensions  and  boundary 
conditions  as  the  E-glass  epoxy  laminate.  The  material 
properties  were  given  in  Equation  5.2  where  the  tensile 
strength  of  epoxy  was  equal  to  54  MPa.  The  laminate  is 
assumed  under  Hertzian  impact  with  projectile  velocity 
equals  40.0  m/second.  Figure  6.5  shows  the  contours  of 
function  f value  for  both  interfaces  at  time  equals  50 
microseconds.  The  f=l  contour  propagations  for  both 
interfaces  are  shown  in  Figure  6.6  for  time  from  50  to  90 
microseconds.  The  same  observations  as  for  E-glass  epoxy 
laminate  are  found  in  the  present  example. 

6.4  Summary  Remarks 

The  delamination  cracks  of  E-glass  epoxy  laminate  and 
graphite  epoxy  laminate  subject  to  localized  impact  have 
been  analyzed  and  presented  graphically.  The  experimental 
result  showed  that  the  fiber  breaking,  transverse  cracks, 
and  delamination  cracks  were  often  observed  in  the  laminated 
plate  when  impacted  by  a projectile.  The  failure  criterion 
only  takes  into  account  the  interlaminar  transverse  shear 
stresses  and  transverse  normal  stress,  which  are  believed  to 
be  responsible  for  the  delamination  crack.  In  the 
computation,  the  laminate  was  assumed  within  the  linear 
elastic  range,  and  the  material  stiffness  did  not  decrease 
after  the  initial  delamination  crack  occurred. 
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(a)  First  interface 


(b)  Second  interface 


Contours  of  failure  value  for  both  interfaces 
of  a 3-layer  graphite  laminate  under  Hertzian 
impact  at  time  equal  to  50  microseconds, 
h=3.81  mm,  v=40  m/second 


Figure  6.5. 
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The  f=l  contour  propagations  for  both  interfaces 
of  a 3-layer  graphite  laminate  under  Hertzian 
impact  for  time  from  50  to  90  microseconds, 
h=3.81  mm,  v=40  m/second 


Figure  6.6. 
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assumption  may  be  one  of  the  reasons  that  the  predicted 
result  could  not  agree  closely  with  the  experimental  result. 
In  the  numerical  results  of  the  three-layer  cross-ply 
laminate  subject  to  impact  loading,  it  was  found  that  the 
delamination  crack  began  near  the  area  of  impact  and  may 
propagate  outwardly  toward  the  edge  of  the  laminate.  The 
results  also  showed  that  the  delamination  crack  for  both 
interfaces  apparently  occurred  at  the  same  time. 

Although  the  results  predicted  herein  could  not  match 
the  experimentally  determined  results  closely  or  describe 
the  actual  phenomenon  of  the  delamination  crack,  they 
provided  some  valuable  information  about  the  initiation  and 
the  characteristics  of  the  delamination  crack  which  agreed 
with  the  observations  of  the  experimental  works.  It  is 
expected  that  this  information  will  be  used  as  a starting 
point  for  the  delamination  crack  studying  and  it  is  hoped 
will  lead  to  further  understanding  of  the  delamination 
propagation  mechanism. 


CHAPTER  7 


SUMMARY  AND  CONCLUDING  REMARKS 


In  view  of  the  increasing  interest  in  using  composite 
materials  for  aerospace  structures,  the  analysis  of 
composite  laminated  plates  becomes  essential.  A systematic 
study  of  the  response  of  laminated  plates  subject  to  static, 
dynamic,  and  impact  loading  by  using  hybrid  stress  finite 
element  method  has  been  conducted.  The  hybrid  stress  model 
developed  herein  is  based  on  the  modified  Reissner 
principle.  The  transverse  shear  effects,  as  well  as  the 
coupling  effects  of  stretching  and  bending,  are  included  in 
the  model.  The  results  of  this  study  showed  that  the  hybrid 


stress  model  is 

sufficient 

and  accurate  to  describe 

the 

laminated 

plate 

response . 

With  the  knowledge 

of 

the 

laminated 

plate 

behavior , 

the  structure  members 

can 

be 

constructed  to  meet  the  specific  design  requirements.  A 
summary  of  present  work  is  described  in  the  following 
paragraphs . 

In  chapter  1,  the  analytical  and  numerical  work  on  the 
response  of  laminated  plates  has  been  reviewed,  with 
emphasis  on  the  hybrid  stress  finite  element  method.  In 
addition,  a brief  introduction  of  laminated  composite 
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structure,  finite  element  method,  and  impact  mechanism  on  a 
laminated  composite  plate  were  summarized. 

In  chapter  2,  the  formulation  of  the  hybrid  stress 
finite  element  method  has  been  described.  The  details  for 
the  derivation  of  the  stiffness  and  mass  matrices  for 
dynamic  problems  were  presented.  The  formulation  for  the 
static  case  were  also  considered.  The  basic  three- 
dimensional  strain-stress  relations  for  an  orthotropic 
lamina  of  arbitrary  orientation  were  examined.  The  terms  of 
the  transformed  reduced  compliance  matrix,  which  were  used 
in  the  finite  element  programming,  have  been  described.  The 
stress  and  displacement  field  of  an  isoparametric 
multilayered  plate  element  were  developed.  All  six  stress 
components  were  included  and  interpolated  within  each 
lamina.  The  stress  field  satisfied  the  homogeneous 
equilibrium  equations  exactly.  The  interlaminar  surface 
traction  continuity  and  upper/lower  laminate  traction  free 
conditions  of  the  multilayered  plate  element  were  enforced. 
All  components  of  the  displacement  field  were  included  and 
interpolated  within  each  lamina  through  shape  functions.  The 
effects  of  transverse  shear  deformation  were  included  by 
allowing  the  rotation  of  cross  section  of  the  individual 
lamina  to  be  independent. 

In  chapter  3,  the  deflection  and  stress  distributions 
for  a laminated  plate  subject  to  static  sinusoidally 
distributed  loading  by  using  the  hybrid  stress  finite 
element  method  were  presented.  The  frontal  solution 
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numerical  technique  was  employed  in  the  solution  of  the 
finite  element  equilibrium  equations.  A comparison  of  the 
results  with  the  available  analytical  solutions  was  made  and 
the  excellent  accuracy  of  the  present  hybrid  stress  model 
was  validated.  Severe  warping  of  cross  sections  was  closely 
predicted  by  the  hybrid  stress  finite  element  method  for  the 
thick  laminate  case.  (i.e.  low  span-to-thickness  ratio). 
Since  the  classical  lamination  theory  ( CLT ) did  not  take 
into  account  the  effects  of  transverse  shear  deformation, 
the  results  predicted  by  the  CLT  were  accurate  only  for  the 
thin  laminate  case.  The  characteristics  of  the  present 
hybrid  stress  model  were  also  discussed.  Fast  convergence 
and  excellent  accuracy  were  observed  in  the  numerical 
results.  For  a thin  laminate,  a more  refined  mesh  was 
required  to  obtain  an  accurate  stress  field. 

In  chapter  4,  the  analysis  of  the  free  vibration  and 
transient  response  of  laminated  plates  by  using  the  hybrid 
stress  finite  element  method  has  been  presented.  A subspace 
iteration  method  was  used  to  compute  the  fundamental 
frequencies  of  a simply  supported  square  laminated  plate. 
The  effect  of  degree  of  orthotropy  of  the  individual  layer 
and  the  effect  of  number  of  layers  on  the  fundamental 
frequency  of  symmetric  and  anti-symmetric  cross-ply 
laminates  were  investigated.  The  results  were  compared  with 
the  3-D  elasticity  solutions  and  excellent  agreement  was 
observed.  The  CLT  overestimated  the  fundamental  frequencies. 
The  effect  of  length-to-thickness  ratio  on  the  fundamental 


164 


frequencies  was  also  examined.  The  transient  response  of  a 
simply  supported  cross-ply  laminate  subject  to  sinusoidally 
distributed  loading  has  been  investigated.  The  Newmark 
direct  integration  method  was  used  to  integrate  the  finite 
element  equations  of  motion  step  by  step.  The  deflections 
and  stresses  as  a function  of  time  were  presented. 

In  chapter  5,  the  dynamic  response  of  laminated  plates 
subject  to  impact  loading  has  been  analyzed.  E-glass  epoxy 
laminates  and  graphite  epoxy  laminates  were  used  in  the 
investigation.  Two  types  of  force  distribution  were  modeled 
in  the  finite  element  programming;  one  is  based  on  the 
principle  of  linear  impulse  and  momentum  with  assumed 
contact  time  and  the  other  is  based  on  the  Hertzian  impact 
law.  The  results  showed  that  the  Hertzian  impact  law  is 
better  in  the  description  of  contact  behavior  and  gives  more 
accurate  solutions.  The  deflection  and  stress  distributions 
along  the  x axis  and  the  y axis  showed  evidence  of  flexural 
wave  motion  induced  by  impact  on  the  laminated  plate.  The 
flexural  wave  propagated  faster  in  the  direction  which  has 
greater  bending  stiffness.  The  interlaminar  normal  and  shear 
stresses  were  dependent  on  the  history  of  the  applied 
loading.  The  deflections  and  stresses  were  proportional  to 
the  projectile  velocity,  but  the  period  of  response  was 
little  affected  by  the  projectile  velocity.  The  effect  of 
laminate  thickness  on  the  transverse  deflection  and  stresses 
has  been  examined.  The  amplitude  and  the  period  of  response 
of  the  central  transverse  deflection  decreased  with 
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increasing  laminate  thickness.  The  flexural  wave  moved 
toward  the  edge  of  the  laminated  plate  faster  in  the  thicker 
laminate  than  in  the  thinner  laminate.  In  addition,  the 
effect  of  laminate  anisotropy  on  the  dynamic  response  of  the 
laminated  plate  subject  to  impact  loading  was  considered. 
The  flexural  wave  propagated  faster  in  the  laminate  which 
has  higher  bending  stiffness. 

In  chapter  6,  the  failure  analysis  of  a laminated  plate 
has  been  investigated.  The  interlaminar  transverse  shear 
stresses  and  transverse  normal  stress,  which  have  been  shown 
to  be  responsible  for  the  delamination  crack  failure,  were 
incorporated  into  a simple  quadratic  failure  criterion  to 
perform  the  analysis  and  examine  the  failure  envelope  at  the 
interfaces.  The  plots  showed  that  the  delamination  crack 
began  near  the  area  of  impact  and  may  propagate  outwardly 
toward  the  edge  of  the  laminate.  The  results  also  showed 
that  the  delamination  crack  for  both  interfaces  appeared  to 
occur  simultaneously.  In  the  present  computation,  the 
stiffness  of  each  lamina  were  assumed  to  be  unchanged  after 
the  initial  delamination  crack  occurred.  This  assumption  may 
cause  the  disagreement  between  the  predicted  results  and  the 
experimental  data.  However,  these  results  provided  some 
valuable  information,  which  could  be  used  as  a starting 
point  for  the  further  understanding  of  the  delamination 
propagation  mechanism. 

Based  upon  the  studies  conducted  thus  far,  the  accuracy 
and  versatility  of  the  hybrid  stress  finite  element  method 
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have  been  ascertained.  It  is  clear  that  the  present  hybrid 
stress  finite  element  model  can  be  easily  implemented  to 
attack  complicated  problems  and  give  more  accurate 
solutions.  The  ultimate  goal  of  this  research  is  to 
investigate  the  dynamic  response  of  laminated  plates  subject 
to  impact  loading  with  emphasis  on  the  delamination  crack 
mechanism.  it  is  hoped  that  the  computer  code  developed 
herein  can  be  employed  in  the  general  applications 
effectively  and  reliably.  Some  works  are  suggested  for 
future  effort  to  make  this  research  more  complete  and 
efficient . 

The  existing  computer  facility  (VAX  11/750)  has  limited 
the  applications  to  laminates  with  a few  layers  only. 
However  in  practice,  the  laminates  are  always  presented  in 
the  form  with  large  number  of  layers.  To  overcome  this 
difficulty,  a global-local  laminate  variational  model  [82] 
can  be  employed  to  reduce  the  computation  time  and  save  the 
computer  core  storage. 

In  the  present  computation,  the  projectile  was  a blunt- 
ended  impactor . The  contact  force  was  governed  by  the 
Hertzian  impact  law  and  assumed  uniformly  distributed  within 
the  impact  area.  A linear  or  higher  order  contact  force 
distribution  around  the  impact  area  is  recommended,  which 
may  achieve  better  descriptions  for  the  stress  and 
displacement  field  of  a laminated  plate  when  impacted  by  a 
blunt-ended  projectile. 
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The  present  model  assumed  that  the  laminate  was  within 
the  linear  elastic  range,  and  that  no  failure  occurs.  In  the 
calculation  for  the  dynamic  response,  the  material  stiffness 
did  not  change  after  the  initial  delamination  crack 
occurred.  A method  that  can  be  used  to  analyze  the 
delamination  crack  propagation  adequately  might  be  to  treat 
the  problem  in  terms  of  fracture  mechanics,  but  this  will 
not  be  easy.  A simple  procedure  might  give  approximate 
results  using  the  following  procedures.  At  each  time  step, 
the  predicted  stresses  are  incorporated  into  a failure 
criterion  to  examine  the  failure  envelope.  If  failure  has 
been  noticed,  remodel  the  finite  element  mesh.  The  failed 
area  is  modeled  by  delamination  crack  elements  and  new 
interpolations  for  the  stress  and  displacement  fields  are 
required  for  those  elements.  Compute  the  stiffness  matrices 
for  the  delamination  crack  elements.  With  the  new  stiffness 

matrices,  the  same  numerical  schemes  as  previously  used  are 

» 

employed  to  solve  the  finite  element  equations  of  motion. 


APPENDIX  A 


STRESS  FIELD  WITH  48  STRESS  PARAMETERS 

°x  = &I3  + [(6i-637)/2-619]x  + 814y  + (62_338)x2/4 

+ [ (33-e39)/2-2e2i]xy  + ( 64-340 )x2y/4  + ?[B22+323X 

+624y+325x2+626xy_3(647'Bll+B4+340)x2y/4J 

a y = 315  + [ (35-341)/2-6lg]y  + B16x  + ( 67-343  )y2/4 

+ t ( P6"e42)/2_2g20]xy  + (68“P44)xy2/4  + ?[S27+328x 

+B29y+33iy2+63oxy_3(e46_3io+38+644)xy2/4i 
az  = ti2{ (1-C )2(2+C ) (B9+B10x+Sliy+B12xy)/4 

+ (1+? )2(2-?) (345+B46x+B47y+B4gxy)/4 
- ( 1-? ) 2 ( 1+? ) ( B 2+B 7+B gX+B 4y ) /4 
+ ( 1-? ) ( 1+? ) 2 ( 338+B43+B44x+B40y )/4 } 
xyz  = t±{(l-?)  (B5+Bgx)/2  + ( 1-?) (1-3?) B?y/8 
+ (1-?) (-1-3?) Bgxy/4  + ( 1 + ?) ( B41+B42x)/2 
+ (1  + ?) (1+3?) B43y/8  + ( 1+?) (-1+3?) B44xy/4 
+ ( 1 - ?2 ) { ( B29+ B3 3 )+ ( B3q+2 B3 3 ) x 
+ [ -3  ( 345-B9  + B2+B38)/4  + ( 33  B2  5 ) ]y 
+ [-3(B46-B10)/2]xy+3(B12-B48)xy2/8}/2} 
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Txz  = ti{ (1-C) (B1+B3y)/2  + (1-C) (1-3C)32x/8 

+ (1-?) (-l-3C)64xy/4  + ( 1+5) (337+B39y)/2 
+ (1  + 5) (1+3C)338x/8  + (1+5) (-l+35)B40xy/4 
+ ( 1-C2) { (B23+634)+(B26+2e36)y 

+ [-3(B45-B9  + B7  + B43)/4+(325-B3x)  ] x 
+ [ ~ 3 ( B47-B11)/2]xy+3( B12"P4  8)x2y/8^/2^ 

Txy  = 617  + e18x  + e19^  + ^20x2  + ^21^  + C {332+333x+e34y 
+B35x2+t-3(B45_B9+B2+338+67+B43)/4-(B25+33i)]xY 

+336y2+3(ei2_648)x2y2/16} 

where  tj.  = (hi  + 1~hi)/2 


APPENDIX  B 


STRESS  FIELD  WITH  55  STRESS  PARAMETERS 

ax  = ®15  + t ^ /2_^21 ^ x + ^16y  + t ^2_^43^/4_^23/2^x2 
+ [ (63-644)/2-2B24]xy  + [ (B5-B46)/4-B26]x2y 

+ [ ( ^4“345 ) /6-B25/3 ] x3  + ? [ B27+B2gX+B29y+B3QX2+B31xy 
-3(B54-B13  + B5  + B46  + 2B10+2B5  3 ) x2y/ 4 ] 

°y  = ^17  + ^6_^47^//2_^20^  + ^18x  + ^ ^8-^4  9 ^ /4_^23/2  ^ 

+ [ (37-648)/2-2B22]xy  + [ ( Bg-B50 )/4-B25 ]xy2 

+ [(B10-B51)/6-B26/3]y3  + ? [ $32+B33X+B34y+B36y2+B35xy 

-3^53_^12  + ^9  + ^50  + 2^4"t’2^45^xy2//,4^ 
az  = ti2{ (1-? )2(2+?) (B11+B12x+B13y+B14xy)/4 

+ ( 1-t-C  ) 2 ( 2-C)  (652+B53x+B54y+B55xy)/4 
- (1-C)2(1+C) [ (B2+B8)+(2B4+B9)x+(B5+2B10)y]/4 
+ (1-C) (l+?)2[ (B43+B49)+(2B45+B50)x+(B46+2B51)y]/4} 

Tyz  = (1-C) (B6+B7x+B1Qy2)/2  + ( l-£ ) ( 1-3C ) B8y/8 

+ (1-C) (-l-3?)B9xy/4  + ( 1+?) (B47+B48x+B51y2)/2 
+ (1+C) (l+3?)B49y/8  + ( 1+?) (-l+3?)B50xy/4 
+ (l-?2) { (B34+B38)+(B35+2B40)x 

+ [ ~ 3 ( ^52-^ll  + ^2  + 343^//4+^36-^30^  ^ Y 
+ [-3(B53-B12)/2-3(B4+B45) ]xy+3 ( B14-B5 5 )xy2/8 }/2 } 
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Txz  = ti{ (1-?) ( B1+63y+34x2)/2  + ( 1-?) ( 1-3? ) B2x/8 

+ ( 1-0  (-1-3?)  B5xy/4  + ( 1 + 0 ( e42+B44Y+B45x2)/2 
+ (1+?) (1+3?)B43x/8  + ( 1 + ?) (-1+3?) B46xy/4 

+ ( 1 — ? ^^328^"339)"^331"*"234l)y 
+ [ "3  ( 352-^ll  + ^8"l"^49^//  4 + ^ 330_33  6 ^ ^ x 
+ [-3(854-613)/2-3(B1q+851) ]xy+3( 4~B55 )x2y/8 } /2 } 

Txy  = 319  + 320x  + 62lY  + 322x2  + 323XY  + 324^  + 325x2Y 
+ 826xy2  + ?{  B37+B38x+B39y+B40x2  + [ -3(  652-33^3  + 62 

+343+38+349)/4_( 330+336) ]xY+34lY2 
+ 3(B14-B55)x2y2/16} 

where  t 3 = (h^  + ^-h.:)/^ 
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